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I . Introductory background topics 



A. The Schrddinger equation 

XXX 

B . Bosons and Fermions 

i. Ground-state wave function of Bosons 

Consider a one dimensional Schrodinger equation for the ground state, which 
may be taken real. We will ask: is this wave function of the same sign 
everywhere? What we mean is: is it all non-negative or all non positive? Of 
course, if it is all non-positive then we can multiply it by -1 and ask again: is it 
all non-negative? Saying it is of a constant sign is the same as saying it has no 
nodes, i.e. it does not go through the x- axis anywhere. 

Exercise: Prove the following theorem: 

Theorem: The ground-state wave function of a ID wave non-degenerate 
Schrodinger equation has no nodes. 

Proof. Let us assume ipipc) is the (real) normalized ground state. Then it is the 
wavefunction minimizes the functional: 



A robust expression for the kinetic energy expectation value 

The kinetic energy is usually defined as T[ip] = —^f™ m ip(x)*\p" \x)dx. Let us 
call this 7\. For continuous and twice differentiable functions that go to zero at 
±00 an equivalent definition, is T[i/j] — — ^^xp' {x)\ 2 dx . Let us call this 



Em = T[ilj]+V[il>] = 



2m 



/_ oo ^'(x) 2 dx + /_ oo ^(x) 2 nx)dx 



(1.1.1) 
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definition T 2 . We now show that T 2 is more robust as it holds for wave 
functions that are discontinuous. To show that 7\ fails, consider the wave 
function ip (x) = e~l*L Notice, it is normalized (ip \ip ) = 1. Let us also 

ft 2 _| I \x\ 

assume for this example that — = 1. Now, for x i/^OO = — e m — (there 
is a discontinuity of VoOO at x = 0) and Vo 00 = V'oOO- Now r x = 
— f- 00 x Po( x T' l Po (x)dx = — 1 which is physically absurd! Kinetic energy must 
be positive! On the other hand, T 2 = /^ o |i/^o(x)| 2 dx = 1. To see that the second 



result is proper, let us smooth the cusp by defining i/j a (x) = e vl*l+«/. Clearly, 
when a — we have are back at the original function \p Q {x), but for positive a, 

no matter how small, the first derivative iO' a (x) = — ^ 2a+ ^ e [ s 

^ av J (|x|+a) 2 

everywhere continuous while xp' (x) has a discontinuity at x = 0. By using 
nonzero a we regularize the discontinuity. Now, we can calculate the integrals 
numerically and we: 



a = 0.1 (i/>a\*l>a) = 1-152 7^] = 0.912 T 2 [^ a ] = 0.912 

a = 0.01 bl> a \il> a ) = 1.019 7^] = 0.988 r 2 [^ a ] = 0.988 
a = 0.001 <tf/ a |tf/ a > = 1.002 7^] = 0.999 r 2 [tf/ a ] = 0.999 
Both T x and T 2 give here the same result now since \p' a is continuous. 

From the table it is clear that lim a ^ 7\[i/; a ] = lim a ^ T 2 [i/j a ] = T 2 [\p Q ~\ 
Tif^o]. Clearly, the expression 7^ is not suited for discontinuities in xp'. 



Using the T 2 expression, we infer that T[ip] = T[\ip\]. Explain why. Is this true 
if ip(x) has a node at x ? Now explain why V[x/j] = V[\ifj\]. 

Now there are two options. Either \xp\ = +xp. In this case there is no sign 
change and there is no node. Otherwise there is a node. Assume it is in x . 
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Let us compute the energy change when one makes a small perturbation to 

p oo 

8E[\xp\]=8T[\xp\]+2 \xp(x)\(V(x)-E[\xp\])S\xp(x)\dx (1.1.2) 

J — oo 

Verify this and show we explicitly assumed that the norm of xp is 1. 
Assume the node is at x = so the wave function changes sign there. 

,g(*) 




Figure 1-1: The absolute value of the wave function and the parabola in it. 

We plot in Figure 1-1 a typical situation. We determine a large enough 
parameter a and the parabola: 

1 

p(x) = -a(x - x' ) 2 + c 

which has the properties that it is tangent at some x x and x 2 to |t/^(x)| (and 
where x 1 < x' < x 2 ). Then one defines a new wave function: 



\ip(x) | x < x 1 or x > x 2 
p(x) x 1 < x < x 2 



(1.1.3) 



By increasing a the parabola becomes narrow; adjusting x' accordingly we 
can cause x x and x 2 to uniformly approach as close as needed. Under these 
conditions, it is possible to show that 
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1^(0)1 



1 |^'(0)| 

c 



2 a 

Now, V[x] — V[ip] oc aT 2 . One power of aT 1 comes from integration over the 
interval of length x 2 — x x where the two functions differ. The second power 
comes from the difference p(x ) — xp(x ) — c oc a" 1 . The kinetic energy 

difference is negative: T[%] —T[\p] « — As a -> oo the total change in 

energy is dominated by the change in the kinetic energy and is thus negative, 
showing that E\j] < E[\xp\] = E[xp]. This is a contradiction, since ip is assumed 
the ground state. ♦ 

The proof can be extended to any number of dimensions. It can be used to 
prove that the ground state of a many-boson wave function has no nodes. 

Prove an important immediate corollary from this theorem: the ground state 
of boson systems is non-degenerate. 

ii. Ground-state wave function of Fermions 

XXX 

C. Why electronic structure is an important 
difficult problem 

In this course, we will study methods to treat electronic structure. This 
problem is often considerably more complicated than nuclear dynamics of 
molecules. The most obvious reason is that usually there are many more 
electrons than nuclei in a molecule. However, there are more reasons: 
electrons in molecules have extreme "quantum mechanical" character, while 
nuclei are more "classical". This last statement requires probably some 
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explanation. We will discuss in this course ways of mapping the electrons in a 
molecule onto a system of non-interacting particles. In such a picture, each 
electron "has" its own orbital or wave-function. Generally the orbitals of 
different electrons in molecules strongly overlap. Furthermore, most of the 
valence orbitals are widely spread over space, encapsulating many molecular 
nuclei. Nuclei on the other hand, are massive, so their wave functions are 
strongly localized and hardly overlap once again, they behave much more 
like classical particles. 

Describing electronic structure is therefore a quantum mechanical problem. 
Now, the electrons interact, via the attractive Coulomb force, not only with 
the stationary nuclei but also with each other: each electron repels each other 
electron via the repulsive Coulomb force. This makes the electronic structure 
problem a many-body problem. 

Quantum many-body problems are very difficult to track down, much more 
difficult than classical mechanical ones. The reason is that for N e electrons 
the many-body wave function is then a function of 3iV e variables, x/j[r v ..r N j, 

where r ; (i=l . . . N e ) is 3 dimensional the position vector. Highly accurate 
numerical representation of such functions is next to impossible for -AT >1. 
Thus, electronic structure methods invariably include approximate methods. 
Typically, the wave function encapsulates much more information than we 



probability for finding an electron at r : and an electron at r 2 and an electron at 
r 3 etc. Indeed, the total probability is 1: 



care to know about. For example, 




the 




(1.3.1) 
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D. The Born-Oppenheimer theory 
i. The adiabatic theorem 

Suppose a Hamiltonian is dependent parameterically on R = (R 1 ,R 2 , ■■■). We 
write this as H[R]. Denote the eigenstates of this Hamiltonian by 

H[R]*Pn[R] = E n [R]ip n [R] (1.4.1) 

Now, suppose we change R with time, so that we have a trajectory, R(t). The 
Hamiltonian will become time dependent: H{t) = H[R(t)]. Suppose the 
system is placed in its ground state at time t = 0, 0(t = 0) = xp o [R(0)]. It will 
evolve according to the TDSE: 

ih<p(t) = H{t)(p{t) (1.4.2) 

Now, the adiabatic theorem says that if R(t) changes very slowly then: 

0(t) = e ">(tty o [fl( t )] (1.4.3) 

Namely, except for a purely TD phase, the evolving state stays the 
instantaneous ground state of the Hamiltian H(t). 

To show this, we use the instantaneous eigenstates to expand 0(t): 

oo 



0(0 = £ a n (t)e-^ EMdT ip n [R(t)] (1.4.4) 

n=0 

We want to plug this expansion into the SE. So let us calculate: 

oo 

n =° oo t d-4-5) 

+ Y j *ne-^ En{T)dT m n [R(t)} 



n=0 

And: 
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tf(t)0 = ^ a n e-y° E ^ dT E n ip n [R(t)] (1.4.6) 

71=0 

Equating the two expressions gives:] 

oo oo 

= £ ind n e-^o^W d ^ n [i?(t)] + ]T a n e-H/o B »W*ifi^[R( t )] (1.4.7) 

n=0 n=0 

Now let us multiply by ^ m and integrate (remember that (ipnlipm) = 8 nm ): 

oo 

= iha m e-^ E ^ dx + £ One-sJo^nWdr^^i^ (1 4 g) 

71=0 

Now, we have: 

= ^WmWn) = hPm^n) + {tmtyn) (1-4-9) 

And since (xpm^n) = WW*/ we find: 

Tmn = (WW = -(^n|W* = ~<m (1.4.10) 

The matrix x mn = (ip m \4>n) is the matrix of "time-dependent non-adiabatic 
couplings (TDNACs)" and it is an anti-Hermitean matrix. Furthermore, for 
n = m 

Re{xp n \xjj n ) = (1.4.11) 

Thus: 

oo 

a m = -TwrnOm ~ ^ e 1 ^^'^' x mn a n (1.4.12) 

We see x mn creates the non-adiabatic transitions, i.e. the transitions out of the 
groud state. 

Let us play a bit with t. Take the time derivative of the TISE: 
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# (O^nW + ff(O0n(O = En(t)lp n (t) + £ n (t)0 n (t) 

Multiply by (i/^ | and get: 



(1.4.13) 



0m # <M + (V"m|Wn) = E n (t)(xp m \xp n ) + E n {xp m \xp n ). (1.4.14) 



Since: (0 rn |tf|0 ri ) = ^(^J^ n ) we have: 



xp m U iPn) = E n (t)(xp m \xp n ) + (E n - E m ){xp m \xjj n ). (1.4.15) 



If n = m then: 



E n {t) = i^ n \fi\^ 



(1.4.16) 



This is called the Hellmann-Feynman theorem, showing that the power is the 
expectation value of the rate of change of the Hamiltonian. If n m: 



= (0m|0n) = 



0m H n 



Eyi E m 



(1.4.17) 



This is called Epstein's theorem, giving an expression for the TDNACs. 
Remember that H[R(t)] depends on t through the positions of the nuclei. 
Thus: 



^-R[R(t)]=R N ^-R[R(t)] 
at oR N 



(1.4.18) 



Thus: 



^ m \^-H[R<it)]\ip n ) n 
L E n [R] — E m [R] RM = Z T ^ A(t) 



(1.4.19) 



where: 



T mn [«]=(0m|^0n) = 



E n [R]-E m [R] 



(1.4.20) 
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We see that the TDNACs depend on the velocity of the nuclei. When nuclei 
are small the couplings out of the ground state are small. The r^ m are called 
the "non-adiabatic couplings (NACs)". We find that as long as E n — E m are not 
zero on the trajectory one can always find small enough R that makes the 
TDNACs as small as we wish. All we need for states n and m to stay 
decoupled is for the following conditions to be met: 

|T mn l = 1^(0^(^(0)1 « 1 (1.4.21) 
Thus we can make the NACs small. So all that is left in Eq. (1.4.12) is d m = 
~ T mm a m- Define the non-dynamical phase as: d nd {t) = i ^T mm {t') dt' (note 6 
is real) and then: 

a m(0 = e Wnd ^ (slow processes — adiabatic process). (1.4.22) 
The total state is 

0(t) = a^COe-stf^to^iMflCt)] = e^CO+edCOty^/Kt)] ( L4 - 23 ) 

1 t 

where 9 d (t) = — - j Q E m (r)dr is called the dynamic phase. It is easy to prove 
that if R(t) traverses a closed loop, the non-dynamical phase depends only on 
that loop and not on the way it is traversed. For example, if we traverse the 
same loop using different velocities, the dynamic phase may change but the 
geometric phase will not. The closed loop geometric angle is called the Berry 
phase. The independence of path is a result of the fact that the non-dynamical 
phase is a line integral, and can be made with no reference to time: 



Electron Density Functional Theory 
© Roi Baer 



Page 15 



Ondit) = i \ (lp m (R[T]Mm(R[T]))dT 

Jo 



dR N 
/ T mm(R) dR N 



4>m[R] )Rn dr 



(1.4.24) 



N 
R(t). 



?(0) 

The dynamical phase for example is not a line integral and its value depends 
not only on the path itself but also on the velocity taken along the way. This 
observation makes the non-dynamical phase a special quantity. Berry has 
shown that this quantity can give us information on the way the Hamiltonian 
is dependent on its parameters. For a real Hamiltonian, for example e l9nd 
around a closed path always equals either 1 or -1. If there is an even number 
of degeneracies enclosed by the path it is 1 and if an odd number - it is -1. 



A path in parameter 
"space" - R 2 



ii. Motivation for the Born-Oppenheimer approximation: classical nuclei 

In a molecule, we can think of the nuclei (having coordinates R) as heavy and 
therefore classical particles which are slowly moving. The electrons are light 
and therefore quantum. They have coordinates r (r includes spin) and they 
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feel the effect of slowly moving nuclei. The electron dynamics is controlled by 
the Hamiltonian: 

H e [R]=f r + V(r,R). (1.4.25) 

The potential V(r, R) describes the Coulombic interactions between electrons 
and nuclei (in atomic units): 

V(r,R) = -) — TT + n / + o / ~~E> ■ (1-4.26) 

*TS\RN-r n \ 2Z r J r ™ 2 ht r nm K ' 

n,N I " 71 1 ni=m N=tM 

In general, one wants to assume that the total energy of the molecule in this 
classical approximation is: 

E = T N + {l> e [R]\8 e [R(t)]U>e[R]) (1-4-27) 

Where xp e [R] is the ground state of the electrons at nuclear configuration R. 
The adiabatic theorem states that this is reasonable as long as nuclei move 
slowly. Thus, the adiabatic theorem allows us to write: 

E = T N + V e (R) (1.4.28) 

where V e (R) is the ground state energy of the electrons at nuclear 
configuration R. This energy is a usual classical mechanical energy and the 
Newton equations of motion apply: 

M N R N = F N = -j^-V e m (1.4.29) 
dR N 

We see, that the adiabatic theorem allows us to consider the nuclei as moving 
in a potential well which is essentially the eigensvalue of the electrons. This in 
essence is the BO approximation when the nuclei are classical. 
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iii. The Born-Oppenheimer approximation in quantum nuclear case 

The classical approach motivates a quantum treatment. We are expecting that 
nuclei will not excite electrons very efficiently. That is the motivation for the 
BO approximation. 

The Born and Oppenheimer development is similar to that of the adiabatic 
theorem, although there are no "external" fields. Suppose we have a system of 
fast particles, with coordinates r = (r 1 ,r 2 , ...) and slow particles, with 
coordinates R = (R lt R 2 , ...)■ The Hamiltonian can be written as: 

H = f R +f r + V(r, R) (1.4.30) 

The Schrodinger equation is: 

Hxp n {r,R)=E n xp n (r,R) (1.4.31) 

Note that we can assume these wave functions are orthogonal: 

i^m) = f f xpn(.r,R)xp m (.r,R)drdR = S nm (1.4.32) 

Now, to proceed, let us first consider the fast r-Hamitonian 

F[R]=f r + V(r,R) (1.4.33) 

In this "fast Hamiltonian" F[R], the slow variables are simply parameters (it 
contains no derivatives with respect to R). Thus, it depends on R 
parametrically. The fast eigenavalue problem is: 

F[R)<p k (r; R) = W k (R)<p k (r; R) (1.4.34) 

The eigenvalues are functions of the parameters R. They are called the 
adiabatic (or Born-Oppenheimer) surfaces. The notation with the semicolon 
between r and R is designed to emphasize that (p k are wave functions in r but 
they depend only parametrically on R. This means, for example that the 
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overlap of the fast eigenstates (called adiabatic states ) involves integration 
only over r, the dynamical coordinates, while R is held fixed, since it is only a 
parameter: 

(0 k O?)|0yOO) = j fe (r; R)<f>j(r; R)dr = S kj (1.4.35) 

We cannot in fact say anything really meaningful about (0 fc (i?)|0y(i?')) when 
R R' (except when R' — R is infinitesimal, but we will not go into this issue 
here). Now, we can expand the "real" wave-function as a linear combination 
of the adiabatic functions: 

*l>n(r,R) = ^(p k {r;R)<S> kn {R) (1.4.36) 

k 

We can do this because for any give R (p k (r;R) span the space of wave 
functions dependent on r. In fact, the expansion coefficients, O n/c (i?) are given 
by: 

fcknOO = f P k (r;R)ip n (r,R)dr (1.4.37) 

Now, let us plug Eq. (1.4.36) into the SE (1.4.31): 

E n f n (r,R) = (f R + P)f n (r,R) 

= ^f R [0 fe (r;i?)O fcn (7?)] +^M4(7?)0 fc (r;7?)cI) fen (i?) (1A38) 

k k 

Note that, since T R = Y.n~ ~ where: P N = —ih- — . We have: 



f R [<p k (.r;R)® kn (R)) 

= O fcn (i?)f fi [0 fc (r;i?)] 

+ Z W ( pN 0fc(r; R) ) ( Pw ° to(i?) ) 



(1.4.39) 



N ' N 



+ <P k (r;R)f R <P kn (R) 
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Multiplying Eq. (1.4.38) by ; (r; R)* and integrating over r gives: 

£(A yft + BjI)O fcn (7?) + (f R + W-(fl)) Oy n (Z?) = E n <t> jn (R) (L4 .40) 



fe 

Where: 



Ay* = | 0y(r;/?)*[f fi fc (r ;j R)]dr = (0y|f R fe ) 
p , _ y j0;(r;/?r[P w fc (r ; fi)]dr _ y (0 y |P w fc ) 

N 

Zh 2 I d \ 3 _ y tv N a 
JhiV J dR~ (pk ldR~: = LW N T i k dR 



N lN ' " N 



The matrices jf k = (<pj ^~4 > kJ are the non-adiabatic couplings in the "fast" 

system. These are exactly the non-adiabatic coefficients in the adiabatic theory 
(Eq. (1.4.20)). 

It is possible to show that: 

£(A, ft + B? k )<t> kn (R) + f R ;n (i?) 



2 , a ,2 



k N 



(1.4.42) 



Thus Eq. (1.4.40) becomes: 

Z Z " 2M~ \W + T *) ° kn(i?) + W jW*JnW = E n®jn(R) (1-4.43) 



k N 



This is a Schrodinger -like equation which determines the coefficients 4> ; - n (i?). 
These are called the "slow eigenfunctions". Once they are computed one has 
an exact solution to the Schrodinger equation. However, we do not really 
want to solve this infinite set of coupled differential equations. Thus we 
assume that the quantities xf k for j k can be neglected. Note that here there 
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is no R N which can be taken as small as needed to make the effect of xf k as 
minute as we need. Still we can jope that the fact that we chose R to be slow 
degrees of freedom allow us to make just this approximation! This was the 
idea of Born and Oppenheimer who neglected xf k : 



The resulting equation is a Schrodinger equation for the slow degrees of 
freedom, which move under the force of electrons derived from a potential 
Wj(jV). When applied to molecules the slow degrees of freedom are usually 
the nuclei and the fast - the electrons. . There is a problem with neglecting rjy 
because of their non-dynamical effect. Taking them into account results in 
treating them as a magnetic field: 



The BO approximation breaks the molecular SE into two consecutive tasks. 
First, the electronic SE Eq. (1.4.34) must be solved for any relevant clamped 
position of nuclei R. Then, the nuclear equation (1.4.26). 

Further reading on this subject: M. Baer, Beyond Born-Oppenheimer: 
electronic non-adiabatic coupling terms and conical intersections (Wiley, 
Hoboken, N.J., 2006). 

E. Electron correlation 

i. The electronic wave function of two non-interacting electrons 

In order to appreciate the complexity of the electronic wave function, let us 
first study a simple system, of two non-interacting electrons in a ID "atomic" 
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(1.4.44) 



Y i 2M~{ i dR~~ Tmm ) ®jnW + Wj(R)4> ]n (R) = E n <t> ]n (R) (1.4.45) 



N 



well. We consider an atomic well given by the potential v ext (x) and we place 
in it an electron. The Hamiltonian is: 

h = -^—^+v ext {x) (1.5.1) 

Now consider a 2-electron problem. Assume we have two electrons, 
Fermions, which are non-interacting in the well. The Hamiltonian is 

H = h(l) + h(2) (1.5.2) 

The notation h{i) means the Hamiltonian of Eq. (1.5.1) with the coordinate of 
the i-th electron. 

What are the eigenstates in this case? First, since each electron can have a 
spin, we must decide on the spin of the state. For now, let us assume the state 
is spin-polarized, i.e. that the total spine is 1, both electrons are in spin-up 
orbitals. We try the following form as a wave function: 

1 

l F(x 1 ,x 2 ) = ^=[i/; (*i)^i(*2) -^o (*2)<M*i)MlM2) (1-5.3) 

Notice that the spatial part is anti-symmetric while the spin part is symmetric. 
This renders the entire wave-function anti-symmetric, in accordance with the 
Pauli principle. The notation a(X)a(T) means both electrons have spin 
projection up (a). We do not yet know if and under what conditions this wave 
function can actually describe eigenstates of the two electrons the 
Hamiltonian (1.5.2). We assume that ip (x) and t/^Cx) are orthonormal. This 
makes the state normalized, since: 
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(1.5.5) 



(1.5.6) 



j H / (x 1 ,x 2 ) 2 dx x dx 2 = i J (^0(^1)^1(^2) - ^o(*2)^i(*i)) 2 dx i dx 2 

= \\ [(^o(xi)^(x 2 )) 2 + (^ (x 2 )^i(x 1 )) 2 ( L5 - 4 ) 
- 2^o(x 1 )^ 1 (x 2 )^ (x 2 )t/; 1 (x 1 )] dx x dx 2 

The first and second terms are both equal to: 

= 1 

The third term in (1.5.4) is: 

j ^o(*l)*M*2)^o(*2)^l(*l) dx l dx 2 

= j ^o(*i)^i(*i) dx i j ^0(^2)^1(^2) dx 2 = x = 

Indeed the entire wave function is orthonormal (thanks to the factor 1/V2 in 
(1.5.3)). 

Now, let us see under what condition the wave function in Eq. (1.5.3) is an 
eigenstate of the Hamiltonian in (1.5.2) 

H¥(x 1( x 2 ) = (fi(l) + h(2))^{x 1 ,x 2 ) 

= -^{i/>i(x 2 )/i(l)i/>o(x 1 ) + 0o(*i)£(2)tfi(* 2 )} (1.5.7) 
- -^{^ o (x 2 )/i(l)0 1 (x 1 ) + t/; 1 (x 1 )/i(2)t/; (x 2 )} 

This should be equated to: 

W(x 1( x 2 ) = -^{£^ (*i)'M*2) - ^0(^2)^1(^1)} (1-5-8) 

If we choose the orbitals ip (. x ) and ^i(x) to be eigenstates of h (which are 
orthogonal so that is compatible with our previous assumption): 
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h\pi{x) = etxpiix), 1 = 0,1,... 



(1.5.9) 



Thus: 



HV>(x 1 ,x 2 ) = ^={e ^ife)^ofe) + ^^0(^1)^1(^2)} 
1 

- ^j={e 1 -ip {x 2 )xp 1 {x 1 ) + e ^i (^1)^0(^2)) 
= (e + ^i) x i , (x 1 ,x 2 ) 



(1.5.10) 



And we see that indeed l F(x 1 ,x 2 ) is an eigenstate with energy 
E = e + 6i. 

ii. Correlation in action: a wave function of 2 interacting particles in an 
Harmonic trap 

We now build a simple electronic structure model that will allow us to study 
in some detail the most basic concepts. For this, we suppose that the electrons 
are in a harmonic atom, that is the potential well: 



(1.6.1) 



The two lowest eigenstates of a Harmonic oscillator are: 



^o(x) = iV e 2 ft 




(1.6.2) 



The normalization constants are: 




(1.6.3) 



And the eigenvalues are: 




(1.6.4) 



Electron Density Functional Theory 
© Roi Baer 



Page 24 



The groundstate energy of the 2-electron system in the triplet state will be 
placing one electron in xp (x) and another in xp^x): 

E = e + e 1 = 2h(0 (1.6.5) 

As discussed above singlet and triplet two-electron ground state wave 
functions composed of two orbitals must be space-symmetric or 
antisymmetric, respectively. We consider below 3 wave functions. The first 
H^ oq is the ground state singlet where both electrons are in ip . The second 
and third are a singlet and a triplet made from one electron in xp and the 
other in xp^. 

1 

*Fs,0l(*l.*2) = ^=[^0(^1)^1(^2) + ^0(^2)^1(^1)] 

= N Ql e-TT^ l+x ^\x 2 + Xl ) (1-6-6) 
1 

^T.Olfrl^z) =^=[^0(^1)^1(^2) -^0(^2)^1(^1)] 

= N 01 e 2h( x i +x 2)(x 2 -x 1 ) 
The normalization factors are: 

N °° = N l = h "» = ^VT^ = iS (1 ' 6 ' 7) 

Eq. (1.6.6) describes the distribution of positions of both electrons in their 
corresponding states. We now want to ask, how much are the electrons in this 
state aware of each other? Do they correlate their motion in some way? One 
way to measure correlation is to consider the quantity (x 1 x 2 ) — (x 1 )(x 2 ). If 
electrons are completely unaware of each other, this quantity, called the 
position-autocorrelation function is zero because then the average of the 
product of their position must decompose to the product of the average. Any 
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deviance from zero indicates some degree of correlation. For the triplet wave 
function in Eq. (1.6.6) we find: 

1 

(*!> = (^r.oi l*i l^oi) = 2 <^o^i - Mo\xi\Mi ~ ^1^0) 
1 

= zWolximMim) + l*i 1^1X^0 l^o) ^ 6g ^ 

1 

= 2«^okil^o) + <^ikil^i)) = 

Of-course, the same result would be obtained if we calculate (x 2 ) because the 
electrons are equivalent. Furthermore: 

. . 1 

<X!X 2 ) = (^r.oiki^ I ^t,oi) = 2<^o^i -^i^ol*i*2l^o^i -^1^0) 

1 

= o «^okl l^oX^l 1*2 l^l) + <^ll*2l^l)^ol*2l^o) 

2 (1.6.9) 

- <^o 1*1 l^lX^l 1*2 l^0> - (^ll^! l^oX^o 1*2 l^l)) 

= -K^ |x|^ 1 )| 2 = -\ — 

2 pLCd 

This negative quantity is there because the Pauli principle pushes the 
electrons to opposite sides (when one electron has positive x coordinate the 
other has negative and vice versa). Let's see what happens in the singlet wave 
function • Here too (x x ) = (x 2 ) = 0. Then: 

{x x x 2 ) = (^00 l*i*2 l^oo) = (^0^01*1*21^0^0) 6 
= <^ol*il^oX^ol*2l^o) = 

Thus, the singlet ground state shows no correlation. However, this does not 
mean that all singlet wave functions of 2 electrons have no correlation. 
Indeed, let us study the situation in W Si01 . The development is very similar to 
Eq. (1.6.9), except that the minus sign is now a plus sign so: 
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<XiX 2 > = ( l F s>0 i|x 1 x 2 | l F s>01 ) = -(xpoxpt + xp 1 xp \x 1 x 2 \xp xp 1 +xp 1 ip ) 



(1.6.11) 



= K^oM^i>l 2 = ~ — 

2 fj.o) 



Here we find positive autocorrelation, indicative of the fact that spin-opposite 
non-interacting pairs of electrons want to stick together (when one has 
positive x the other wants to have as well and when one has negative the 
other one wants negative as well) like non-interacting bosons. 

Since there is no interaction between the electrons, the correlation in these 
wave functions arises only from the Pauli principle, i.e. because we impose 
the fact that electrons are Fermions. This is called Fermi correlation. Our 
lesson is this: 

1) Wavefunctions that are mere products of singe-particle orbitals have no correlation. 

2) If the products are symmetrized like in the case of the excited singlet the correlation 
{x 1 x 2 ) is positive indicating that the particles "like to be together" i.e. both on the 
right or both on the left of the origin. 

3) If the products are anti-symmetrized like in the case of the triplet the correlation 
{x 1 x 2 ) is negative indicating that the particles "like to be seperated". 

Up to now, we assumed no e-e interaction. So now let's include it and add to 
the Hamiltonian an interaction term: 



Let us take one case of coupling which is simple enough to yield to analytical 
analysis: 



With y 2 < co 2 . This interaction seems strange at first site because it does not 
depend on the distance between the particles, as we are used to from 
electrostatics, yet, it does describe a repulsion: since if x x and x 2 are both large 
and of the same sign this is energy-costly; if they are both large and of 
opposite sign that lowers energy. In this case, the Hamiltonian is: 



H = h(X) + h(2) + V(x lt x 2 ) 



(1.6.12) 



V(x x ,x 2 ) = ^7 2 x 1 x 2 



(1.6.13) 
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(1.6.14) 



Q Pi . P2 1 22,1 92, 2 

Exercise 1-1 

Find the eigenvalues and eigenfunctions of this Hamiltonian. 

X-hx X—x 

Solution: Define new coordinates X and x by: -j=- = x x and = x 2 . The 
conjugate momenta are: p x = -j=- and p 2 = — j=- (show that this is indeed so by 

calculating the commutation relations [P.A'], [P, x] etc.). Then the new 
Hamiltonian is: 



^ (P+p) 2 (P-p) 2 1 7 fX + x\ 2 1 7 (X-x 



//= ~ ~~ + ~4iJ~ + 2 fia) \ V2 y '2-VV2 



4/i 



/A + X\" 1 (A — XV 

br) + 2^ 2 hr) 



_ /X + x\ fX — x\ 

^ (—)(—) 



(1.6.15) 



Or, after rearranging: 


H = 




+ 


p 2 1 

^ + 2^ 2 -K 2 )x 2 




(1.6.16) 



We see that X and x do not interact and each is a Harmonic oscillator on its 
own. Let us define cos 6 = and H x = ^/co 2 + y 2 = V2co cos ~ and H 2 = 
-y/o) 2 — y 2 = V2o> sin ^. We find: 



W nm {X,x) =\f; n (X;Cl 1 ,ii)\lt Tn (x;Cl 2 ,ii), 



(1.6.17) 



With xfj n (z; £1, fi) = N n e 2H n _ 1 (q) / n = 0,1, q = l—z, H n (q) are Hermite 



polynomials of order n and iV n are orthonormalization constants. 
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Exercise 1-2 

Write the ground-state energy and wave function for the triplet state of the 
system in the previous exercise. Determine the effect of interaction on the 
energy by calculating r = — and on the correlation function c = (x 1 x 2 ) — 

Egs,y=0 

{x 1 ){x 2 ). 

Solution: We need to impose upon the spatial wave function of (1.6.17) to be 
antisymmetric. The two variables X = Xl ^ 2 and x = Xl ^ 2 are respectively 

symmetric and antisymmetric combinations of the positions of the electrons. 
Since x is the antisymmetric combination we require the Hermite polynomial 
in x to be odd. The lowest antisymmetric state is the combination n = for X 
and m = 1 for x: 



HfliX 2 iiCl 2 x 2 



Wo^X.x) = N 01 e 2h 2ft x 

1 3 1 ,_/ 9 G\ 

£"„, = —Mli + - Ml, = -/io) x V2 cos - + 3 sin- 
gs 2 1 2 2 2 V 2 2/ 

Now, let us write the wave function of Eq. (1.6.18) in terms of x x and x 2 : 



(1.6.18) 



_^kCr +r V ^ 2 (r -r V ( X l — X 2) 

V2 

= -£±e is - (*i + *2>e 2ft— '^(x, - x 2 ) 
V2 



(1.6.19) 



One can compare the effect of the interaction by looking at the ratio between 
ground-state energies of the system with and without interaction as a function 
of 6, the interaction strength: 

E -fl 1 +-fl 3 cos- + 3sin- 

r _ c gs, Y _ 2 1 2 2 _ 2 2 (1.6.20) 

E gs 2o) 2V2 

The result is shown in Figure 1-2 (left panel). 
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n n 

Figure 1-2: (Left panel) The ratio r = E 12i y/E 12 ; (Right panel) the correlation ratio c(G)/ 
c f^J vs interaction strengths cos 9 = J 

The interaction lowers the energy, because now the wave function can acquire 

a structure that promotes the electrons being away from each other. Thus one 

is pushed towards the +x direction and the other towards that of — x and thus 

they acquire a large negative value of x 1 x 2 . To see this note that the 

expectation values of x and X are both zero and therefore (x^) and (x 2 ) are 

zero as well. Furthermore note that X 2 — x 2 = 2x x x 2 and the auto-correlation 

is: 

c(0) = ( Xl x 2 ) - ( Xl )(x 2 ) = ( Xl x 2 ) = \{{X 2 ) - (x 2 )) (1.6.21) 

The expectation value of the square position in harmonic oscillator is easily 
obtained using Hellman-Feynman theorem (rederive Eq. (1.4.16) in terms of 
any parameter-dependent Hamiltonian, not necessarily the time t): 

H = T + ^a> 2 y 2 (n + ^h = d ^ = fico(y 2 ) n (1.6.22) 

-ft -ft 
So: (X 2 ) = and similarly {x 2 ) 1 = Thus: 
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c(£0 = ( Xl x 2 ) - ( Xl )(x 2 ) = ( XlX2 ) = \(X*- x*) = i£ gi- - § 



(1.6.23) 



e . e 
cos- sin- 



For no interaction (6 = -), there is only correlation due to the Pauli principle: 
since both electrons have spin up they cannot occupy the same point in space) 
We see, as derived above: 

/7T\ h 

c I — ) = Fermi Correlation = (1.6.24) 

The ratio between full correlations and Fermi correlation is shown in Figure 
1-2 (right panel). The correlation is negative and larger (in absolute value) than 
the mere Pauli correlation. Indeed, the interaction pushes electrons away from 
each other. 

Looking at the wave function of Eq. (1.6.19), it is evident that because n x fi 2 
there is no way to write this as an antisymmetrized sum of products of 1- 
electon functions. From these exercises with harmonic oscillator systems we 
find that the issue of correlation can be quiet complicated. In realistic 
electronic systems, when the interaction is Coulombic, not Harmonic, the 
situation is even more complicated because of the lack of analytical solutions. 

F . The electron density is a much simpler object 
than the wave function 

The complexity of the wave function is overwhelming. It includes all the 

information we can have on the molecule in a certain state. However, all these 

intricacies and details are often uninteresting for us: in many cases, we simply 

have no direct use for them. Take, for instance, the electronic energy - our 

primary target in the Born-Oppenheimer picture. It depends only on the 
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relative distance of pairs of particles. This is because the electron-electron 
Coulomb repulsion is a pair-wise interaction. 

One interesting quantity, besides energy, that can be extracted from the 
electronic wave function is the electron density „ Cr 5 . This 3D function tells us 
the expectation value of the density of electrons. That is n(r)d r is the 

expectation value of the number of electrons at a small volume d 3 r around 
point r . Thus, we can write: 

n(r) = (xp\n(r)\xp) (1.7.1) 

We use the notation: 

(ip\(j)) = j ip*(r 1 ,r 2 ,...,r Ne )(p(r 1 ,r 2 ,...,r Ne )d 3 r 1 d 3 r 2 ...d 3 r Ne (1.7.2) 

Here n(r) is the operator corresponding to the electron number density. Since 
electrons are point particles, and the position operator for the i th electron is r t 
this operator is defined by: 



n( r ) = Y S(fi - r) (1.7.3) 



i=l 



We used the definition of a 5-function, according to which: 

j S(r 1 -r)f(r 1 )d 3 r 1 =fCr) (1.7.4) 

The "function" 8(r j — r) is the density of electron i at r. 
Exercise 1-3 

Calculate the expectation value of the electron number density and show that: 

n(r) = N e J xp(r,r 2 ,...,r Ne )xp(r,r 2 ,...,r Ne )d 3 r 2 ...d 3 r Ne (1.7.5) 
Solution: Plugging (1.7.3) into (1.7.1) gives: 
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n 



(r) = | 0(: 



r Ne )d i r 1 d 2 r 2 ...d 3 r Ne 




xp(r 1 ,r 2 



,...,r Ne )8(fi 



-r)^(r 1( r 2 ,..., 



(1.7.6) 




-.»w e )S(ri 



r)^(r 1( r 2 ,..., 



The last step stems from the Pauli principle: all electrons are identical, so we 
can replace electron r £ by electron r x . The sum is now over identical numbers 
so it becomes a mere multiplication as in Eq. (1.7.5). 

Looking at Eq. (1.7.5), we see that n(r) involves integrating out a huge 
amount of wave-function details. It is as if only the data concerning the 
density distribution of a single electron remains! This is multiplied by N e in 
Eq. (1.7.5) so n(r) accounts for the combined density of all electrons. Indeed, 
integrating over the entire space, one obtains from (1.7.5): 



Expressing the fact that the total number of electrons is N e . 
Exercise 1-4 

Calculate the ID electron density of the triplet ground state from Exercise 1-2. 
Solution: If the wavefunction of Eq. (1.6.6) is taken then: 




(1.7.7) 
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/-°° ^Ol(x 2 +X2) 

n(x) = N I e a (x 2 — x) 2 dx 

J — CO 



= Ne fi~ 



£ 



2 



e ft xfdx 2 +x 2 J e ft dx 2 



[ 



2 

\lCiiX2 



= N 



— e ft dx 2 [(x^) + x z ] 

m J -co 



= Ne ft 



x 2 \hit\ h 




We choose iV to ensure that j n(r)d 3 r = 2: 



Defining the "average frequency": 

- — ) = 



1 _ 1/ 1 
n ~ 2 vai 





1 cos- + sin- 

-""2 2 



V2co sin 6 
We find the density of the state in Eq. (1.6.19) is: 



n(x) = JV I 

J— c 



fll(x+x2) 2 n2( x_;c 2) 2 

e 2 2 (x — x 2 ) 2 dx 2 



a 3 

!7rH 2 



e - n *'(i + 2(2n 1 -n)x 2 ) 



(1.7.8) 




(1.7.9) 



(1.7.10) 



(1.7.11) 



G. The variational principle 

When we look for the ground-state energy of a complicated system, with 
Hamiltonian H, the variational principle is often extremely important. It says, 
quite simply that the ground-state energy is a minimum of a functional of 
wave functions. The functional is: 



E[<p] 



(0lg|0) 

(010) 



(1.8.1) 
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And the variational principle states that: 

E g s = E[4> gs ]<E[4>} (1.8.2) 

For all 0. What this principle allows us is to determine which of the two wave 
functions X and 2 gives an energy closer to E gs , without knowing E gs . 
Simple the one with lower energy is closer. 

We can build a parameterized family of functions y and determine the 
"optimized" parameter as the one that gives minimum to E(y) = E[(p Y ]. 

Exercise 1-5 

Consider the quantum ground-state problem for a particle in a well v{pc) = 
^/cx 4 . We consider the family of functions: 

a (x) = -= (1-8-3) 

V V27TO" 



Here a is a positive parameter. What is the best function for representing the 
ground state? 

Solution: The functions are normalized. The energy is H — T + V. We have: 



E{a) = E[4> a ] = (0 ff |f |0 ff ) + (4> a V(f> a ) = + \ka* 




(1.8.4) 


Where we used the fact: 


h 2 f°° h 2 




(1.8.5) 


Thus the minimum energy is obtained from: 


h 2 ( h 2 \ 1/6 




(1.8.6) 



And so: 
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2/3 

rr :j i ir 

E* = 

h 2 vl/3 



8m (J*-) 

\12kmJ 



3 / h 2 Y ,D _ 3h 2 l 

+ 4 k \12km) ~16m"^ ( L8 - 7 ) 



H. Dilation relations 

Consider the normalized xp(x). From it, we can define a family of normalized 
wavefunctions: 

ip Y (x) = Y 1/2 ip(.Y*) (1-9.1) 
We can check that each member is indeed normalized: 

(VvhM = / kO)| 2 d* = r J l^(rx)| 2 dx = | |tf/(y)| 2 dy (1.9.2) 

This operation of "stretching" the argument of a function is called dilation. 
Dilation affects functionals. Consider the kinetic energy functional: 

T[xp] = {xp\f\xp) = j ip(x) f-i A_U( x )dx (1.9.3) 

Exercise 1-6 

Prove the dilation relation of T[ip Y ] and T[ip]: 

Solution: We show that this is a simple relation: 
T[xP Y ] = j xP Y (x)(-^-^\iJJ Y (x)dx 

= Y j i{j(yx) [~\^\ *l>(yx)dx (1.9.5) 



= y 2 j $(y)(-\-j-i\M>(y)dy =Y 2 Tbp] 



The potential between all particles is called "homogeneous of order n" if: 

V(yx! yx N )=y n V{x 1 x N ) (1.9.6) 
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Examples are n = 2 is for Harmonic potential wells and harmonic interactions 
and n = — 1 is for the Coulomb potential wells and Coulombic interactions. 

Exercise 1-7 

Prove the dilation relation for homogeneous potentials of order n: 



V[ip Y ] = Y - n V[ip] (1.9.7) 

Solution: For such a potential, 

f 2 
V\$y\ = I Y N \^ Y (y^i'-'Y^N)\ V{x 1 ,...,x N )dx 1 ...dx N = 

= j |^y(yi.-.yiv)| 2 ^(y, ->y)dy 1 ...dy h 



Y Y> (1.9.8) 

= Y~ n j \4> Y (yi- -.yw)| 2 ^(yi--.yw)dyi -dy N 
= r" n W 1 ]- 

We combine the results from Eqs. (1.9.4) and (1.9.7) and obtain an interesting 
property of the total energy for systems with homogeneous interactions: 

E[xp Y ] = T[4, Y ] + Vtyy] = Y 2 T[ip] + Y~ n Vm (1-9.9) 

For a molecule, the interaction between the particles (electrons and nuclei) is 
the Coulomb potential V(r;R) = V(r 1 ,...,r Ne ;R 1 ,...,R NN ) which is 
homogeneous of order n = — 1 one finds the energy of a molecule obeys: 

E[xp Y ] = T[xP Y ] + V[xP Y ] = y 2 7^] + (Coulomb) (1.9.10) 

i. The concept of the virial in classical mechanics 

First, let us define the virial. For a system with coordinate q n collectively 
denoted as q the virial in classical mechanics is the time average of q • F 
where F is the force vector: 

1 f t+ l 

virial = (q ■ F) = - q(t') • F(t')dt' 
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It can be shown, that for bound systems: [5] 

virial = (q ■ F) = —2(T) 

For conservative systems the force is a gradient of the potential, F = —Vv(q). 
The viral relates to dilation of the coordinates through: 

d 



dy 



v Y (q) = q ■ Vv(yq) = q • Vv y (q). (1.9.11) 



For homogeneous potentials we have: y — v y (q) = ny n v(q) = nv y (q), thus: 

dy 

rtVy(q) = yq • Vv y (q) (v homogeneous order n) (1.9.12) 
In particular, for y = 1: 

nv(q) = q ■ Vv(q) (v homogeneous order n) (1.9.13) 
We will especially be interested in Coulomb systems, where n = — 1, then: 

v(q) + q ■ Vu(qf) = (v homogeneous order n = —1) (1.9.14) 
Exercise 1-8 

Eq. (1.9.13) is known as Euler's equation, after its discoverer. In 
Thermodynamics it is extremely useful. Thermodynamics stipulates that if 
you know the energy of a system t/(S, V, iV) as a function of its macroscopic 
parameters S, V, N then you have complete thermodynamics knowledge. The 
first and second laws of thermodynamics state that: 

dU = TdS -pdV + fl- dN, (1.9.15) 

where T, p and jl are respectively the temperature, the pressure and the 
chemical potentials. A second stipulation is that U is a homogeneous function 
of order 1. 

From this, show that for all thermodynamical systems 

1. U = TS-PV + p.-N 
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For homogeneous potentials the viral theorem becomes: 

virial = -n(v) = -2(T) (1.9.16) 
For Coulomb systems n = — 1: 

(v) = -2(T) (1.9.17) 

From which: 

E = (v) + (T) = — (T) = ^ (1.9.18) 
We now show that this relation holds also in quantum mechanics, 
ii. The Virial Theorem in quantum mechanics 

Now, if ip is the ground-state energy then £"[i/; y=1 ] obtains the minimum and 
therefore: 

^E[4> Y ] = (1.9.19) 

Plugging this into (1.9.10), one obtains: 

= 2yT[\l)] - ny-^M^] (1.9.20) 

Since we know that the optimal value of y is 1, then (dropping the [ip] symbol 

and remembering that all following quantities are ground state expectation 
value: 

2T = nV 

2^ m x (1.9.21) 



£=7+7=1+ 



These relations show that the virial theorem (Eq. (1.9.18) holds in quantum 
mechanics provided xp is the full molecular eigenstate. For molecules n = — 1: 
and one has: 
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2T = -V (Coulomb) 



1 

T + V = -T = -V 



(1.9.22) 



E = 



A subtle issue: all energies and potentials in the above expressions are 
absolute. We usually give energy and potentials only to within an additive 
constant. However, the fact that the potential is homogenous, it cannot 
tolerate addition of a constant (x 2 is homogeneous of order 2 but x 2 + a is 
not). 

iii. Slater's theory of the chemical bond using the viral theorem 

What happens when xp is the electronic eigenstate in the Born Oppenheimer 
approximation? If we look only at the electronic wave function we do not 
expect Eq. (1.9.21) to be valid. Indeed, using H e = f + V(r; R) where T = 
( — - V 2 . J and V(r; R) = (with obvious summation on nuclei and electrons) 



Upon dilation x(j Y [R](r) = y 3We / 2 i/;[J?](yr) and note that R is not dilated. Then 
we have for the KE and potential: 




H e [RMR](r) =E[R]ip[R](f) 



(1.9.23) 




2 T[R] 



(1.9.24) 



And we define: 
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l3N er 



w y [r,r'] = (ip Y [R]\v[R']\M R V = f M r > R ¥ v (r;W 

= Y 3N e f xp(y r] Ryv(r,R')d 3N er 

(1.9.25) 

= 7 J \p(s,R) 2 V(s;YR')d 3N es = yW[R,yR'] 

We have defined W, a two J? quantity. Of course, the physical interaction 
potential is V[R] = W[R, R]. The reason we defined W this way is that we can 
now write: 

E Y [R] = Y 2 T e [R] + YW[R,yR] (1-9.26) 

Note that we inserted the subscript e to the KE since now it is only the KE of 
the electrons. 

One still has the variational principle, i.e. dE Y [R]/dY\ Y * =1 = (Eq. (1.9.19)): 

= 2 Y T e [R] + W[R, Y R] + y(R ■ V r ,W)[R,R'] r , =yR (1.9.27) 
Here we used the fact that ^-f(yR) = [R ■ Vf](yR). Putting in the 

u,y 

solution y* = 1 we find: 

= 2T e [R] + V[R] + R • V R >W[R,R'] R=R , (1.9.28) 

From the Hellmann-Feynman theorem (see section XXX) V R '(W[R, R']) R r= R = 
V R E(R). Adding to the electronic energy E the nuclear potential energy 

E N = -Yii*v i„ ' „ i gives the BO energy E B0 = E + E N and since the nuclear 
energy is homogeneous of order -1, we can use Eq. (1.9.14) and: 

= E N + 2^ R / ' Vrj En> (1.9.29) 

; 

Plugging all this into eq. (1.9.28) we find Slater's formula[6]: 
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T e = -E B0 ~ R • VrE B0 (1-9.30) 

From this, using E B0 (J?) = T e + V tot (where V tot = V + E N is the total PE of the 
molecule) we find also: 

V tot = 2E BO + R-V R E BO (1.9.31) 

We point out an interesting property of Eq. (1.9.30) is that we can obtain 
electronic properties such as kinetic energy T e or potential energy V directly 
from the potential surface E B0 . This will be important for analyzing properties 
of DFT quantities. Note that at stationary points on the Born Oppenheimer 
PES where V R E BO =0 the usual virial relation eq. (1.9.21) holds (with neglect of 
nuclear kinetic energies). 

Slater derived Eq. (1.9.30) in a different manner, not from variational or 
dilation considerations. He used it for analyzing the general nature of 
chemical bonding of a diatomic molecule. We follow in his footspteps here. 
Consider a diatomic molecule. Suppose the BO curve is give as a potential of 
the following form, in terms of the inter-nuclear distance R: 

Where A and B are positive constants and typically m » n.For example, in the 
case the PE curve is given by the Lennard-Jones potential then m — 12 and 
n = 6. This potential describes long-range attraction and short-range 
repulsion. It has a minimum at: 

One then has: 
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fl(l-n) _ >t(l-m) _ 

When R -» oo we have: T e « — . Let us consider what happens to 

the KE as the distance R is reduced. Clearly it initially descends, however, at 
some distance R T it reaches a minimum, and below this distance it begins to 
increase. Since: 

d Bnjl - n) Amjl - m) (1 9 35) 

di? e R n+1 R m+1 

We find the minimum kinetic energy is obtained at a distance somewhat 
larger than R Min : 



_( Am(m-l) y-n_,m- ly-n 
Rt ~ Un(n - 1) ) "U-lJ ^->«Min 



For example, for m = 12, n = 6 we find: R T = 1.14R Min . The sum of electronic 
and nuclear potential energy is: 

4(2 -m) B(2-n) „ n 

7 = — i - -1 L + 2E m (1.9.37) 

At large R, V + E N actually increases as the atoms approach. The potential 
energy does not rise indefinitely. At some inter-nuclear distance, larger than 
the distance at which the BO potential obtains its minimum, it reaches a 
maximum and then starts decreasing. We have: 

Am(2 — m) Bn(2 — n) 

^_ ^_ (1-9.38) 

/Am(m — 2)\ m - n /m — 2\m- n 



Rv ~ { Bn{n - 2) ) " U^l) Rmn > Rmn 

For the Lennard Jones potential we have: R v = 1.16R Min . For R < R T we find 
that the kinetic energy grows sharply as R decreases, while the potential 
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energy decreases. This sharp behavior can be interpreted as a condensation of 
the electron cloud around each nuclei as such a condensation causes increase 
in kinetic energy T and decrease in electronic potential energy V. The total BO 
energy continues to drop as long as R > R M in an d then starts rising. 

As an application of this theory consider the X ... Y + system where X is a 
neutral atom and Y + a distant atomic ion (can also be a molecular ion). The 

cc(X} 

energy at large distance R is that of Eq. (1.9.32) with A = and B = where 

<x{X) is the polarizability of X and n = 4. This form of the potential surface 
results from a simple exercise in electrostatics. Quantum mechanics is needed 
only for calculating the value of the polarizability of X, a(X). Slater's analysis 
states that: 

3a(X) 

T = — — — E 

2R 4 

(1.9.39) 

_ a W 
v tot - R4 + 2£"oo 

So the force due to the total potential energy is repulsive while that due to the 
KE is attractive. The KE attraction outweighs the PE repulsion by a factor of 
three halves. Slater's analysis shows that the interaction here is fundamentally 
different from that of electrostatics, where the force attributed entirely to the 
Coulombic attraction. 

Another application is to the simplest chemical bond, that of the molecule. 
The text-book analysis of this molecule assumes that the molecular orbital is a 
"bonding" linear combination of the two atomic orbitals (LCAO), both having 

r 

the Is shape: ip ls (r) = e a o, but each centered on its corresponding nucleus. 
Quantitatively this picture is very inaccurate. While it does show binding, its 
predicted a bond length is almost 1.4 A (while the true bond length is 1.06A) 
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and its estimate of atomization energy is 1.7 eV while the experimental value 
is close to 2.7eV. There are no electron correlation effects here, so the only 
source of inaccuracy is the limited efficacy of the small basis. 




Figure 1-3: The optimal dilation parameter % t as a function of the inter-nuclear distance R 
in H£ , where the wave function is the sum of 2 ls-type orbitals (localized on each of the 
two nuclei). 

One can add a dilation factor to the wave function i/^ s (r) = e a o and perform 
a variational determination of <\ Then, one can obtain, for each value of the 
inter-nuclear distance the optimal denoted At R -* oo the dilation 

parameter should be equal to 1, as then the Is orbital is exact. At R -> the 
molecule becomes the He ion and <"* should approach the exact value of 2. See 
Figure 1-3. But the optimal value <"* is not monotonic! At large R it becomes 
smaller than 1, indicating that the electron cloud expands, the kinetic energy 
drops and the potential energy rises. Thus long-range attraction is driven by 
the kinetic energy. At R = 5.2a the value of <"* is again 1 and then grows 
monotonically as R is reduced. Now the orbital is contracting, kinetic energy is 
growing and potential energy is decreasing. In this range, chemical binding is 
driven by the potential energies. At the distance minimizing the BO energy 
the value of ( is: <"* « 1.23, i.e. the orbitals are then considerably contracted 
(relative to the free H atom). With such contracted orbitals, the calculated 
bond length of the LCAO calculation is 1.06A. The bond energy, computed by 
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subtracting the calculated H2 + energy from that of the hydrogen atom gives the 
bond energy of 2.35eV. The first to report this type of calculation is in ref [7], a 
clear exposition of the calculation is given in refs [8]. These improved 
theoretical estimates, with a two-function basis, achieved with just a single 
dilation parameter, show that orbital shapes are crucial for a good description 
of chemical bonding. 

This result contradicts our intuition regarding the behavior of electrons in a 
covalent bond. It seems that a major source of the chemical bond energy in 
is not due to the electrostatic benefit resulting from the placement of electrons 
in the middle region separating the two nuclei. Rather, it is due to the 
contraction of the electron wave function near each of the nuclei induced by 
change in Coulomb potential near that nucleus due to the presence of the 
other nuclei. Around one nucleus of H2 the nuclear Coulombic potential 
becomes deeper and distorted towards the second nucleus. As a result, a 
contraction of the electron density around the nucleus is obtained, and the 
charge placed between the atoms is not in the "midway" region of the bond: it 
is between them alright but very close to each nucleus. Each nuclei seems to 
share the electron not by placing it midway between them, but rather by 
having the electrons much closer to each nucleus. 

Combining this with the above result, we find a very general conclusion: 

Rule 1 of chemical attraction : As atoms approach from afar, (for R > R T ) the 
kinetic energy decreases and the potential energy increases and thus it is the 
lowering of kinetic energy which is responsible for the distant attraction between 
molecules. 

Rule 2 of chemical bonding : For distances near but larger than the bonding 
distance the attraction is more subtle. The electronic kinetic energy rises 
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sharply and the potential energy drops even more sharply. Thus, a major 
source of energy is the contraction of orbitals around their respective nuclei, 
inspired by the deepening of potential wells due to the presence of close-by 
nuclei. The electronic PE is thus the major glue when atoms are close to each 
other: It has to offset the nuclear repulsion and the repulsion due to kinetic 
energy of electrons. 

iv. Other Dilation Facts 

Consider the Schrodinger equation: 
h 2 

- — ip"(x) + v(x)ip(x) = Eip(x) (1.9.40) 
2m 

We now dilate: we assume ifixOc) — ip(Xx) is an eigenstate and find the 
corresponding Hamiltonian. Note that ^'(x) = X 2 xp"{Xx) thus: 

- ^ 00 + v&xWxto = Exp A (x) (1-9.41) 

We see that the required Hamiltonian involves a dilated potential and a scalled 
mass by a factor A 2 .The energy is left intact. Alternatively we can multiply 
through by the square of the dilation factor A 2 and obtain: 

h 2 

+ X 2 v(Xx)ip x (x) = X 2 Ei/j a (x) (1-9.42) 

2m 

Now the potential is dilated and scaled: v(x) -» X 2 v{Xx) when then the wave 
function is dilated by ip(x) -* xp(Xx) and the total energy is also scaled by 
E -> X 2 E. This is general for basically any Schrodinger equations and holds for 
any number of spatial dimensions. 

Now consider homogeneous potentials in 3D: u(lr) = X n v{r). In this case the 
potential transform is indeed a multiplicative scaling: v(r) -* X 2+n v(r). In 
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particular, in Coulomb systems, when n = — 1, we find that the potential 
scaling is simply by a factor X. We can turn around the discussion and state: 



If the Coulomb coupling constant is scaled by a factor X then the total 
energy is scaled by a factor X 2 and the wave function is dilated by a 
factor X. In particular, the density is dilated and scaled as n(r) -> 
X 3 n(Xr). 

According to the virial theorem, the potential energy and the total energy are 
related by ~ = — T = E, thus the expectation value of the potential and kinetic 
energies is scaled by a factor X 2 as well. 
Exercise 1-9 

An application of the above discussion allows an interesting exact conclusion 
concerning the ground-state energy per particle of the homogeneous electron 
gas of density n (see section XXX for details). Denote this quantity, e HEG (n, e 2 ) 
and the Coulomb coupling constant e 2 /4ne . Show that: 

3n^— + e 2 ^— - = 2e HEG (1-9.43) 
on ae l 

Hint: Start from e HEG (X 3 n, Xe 2 ) = X 2 € HEG (n, e 2 ). Then take the derivative with 
respect to X. Then set X = 1. 

I. Introduction to Functionals 
i. Definition of a functional and some examples 

A formula that accepts a function and maps it to a number is called a 
functional. We denote the number a functional F maps to a function n(r) by 
this by the symbol F[n]. Take for example the simplest functional: F[n] = 0. 
This functional maps any function n(r) to the number zero. A more 
interesting functional is: 

Electron Density Functional Theory Page 48 

© Roi Baer 



F ro [n\ =n(r ) (1.10.1) 

This functional maps to each function n(r) its value at a certain point r . 
Next, consider the functional 

WM = V A n(r)d 3 r, (1.10.2) 

mapping each function its average value in a given volume V. 

Another familiar example is the kinetic energy functional for a particle of 
mass pL in some wave function xp(r): 

h 2 r 

T[ifj] = - — ^*(r)V 2 ^(r)d 3 r, (1.10.3) 

Jy 

Or the energy for a given potential V(r): 

E[xp]=T[xp]+V[xp], (1.10.4) 

Where: 

V[xp] = j V(r)\xp(r)\ 2 d 3 r, (1.10.5) 
ii. Linear functionals 

A functional is called linear if for any two functions n(r) and m(r) one has: 
F[n(r) + m(r)] = F[n(r)] + F[m(r)] and if F[an(r)] = aF[n(r)] for any 
number a. The functionals above /VoW an< ^ hi n ] are examples of such 
functionals, please check. Any linear functional can be written in the form: 

F[n] = j f L (r)n(r)d 3 r. (1.10.6) 

The function f L (r) is called the linear response kernel of F. For example, the 
linear response kernel of F rQ is f L (r) — 8(r — r ). And that of I v is f L (r) = -. 
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iii. Functional Derivatives 

Let us now introduce the concept of a "functional derivative". Consider the 

SF 

functional F ro [n] of Eq. (1.10.1). A functional derivative, denoted 

measures how F[n] changes when we make a small localized change of n(r) 
at F. 

Let us give a more detailed definition. Suppose we change the argument 
function in F[n] by making an arbitrarily small change in (r) -> n(r) + eO(r). 
We use an arbitrary function O(r) multiplied by an infinitesimal number e 
(that is, a number e we intend to make vanishingly small at a later statge). In 

SF 

principle, for any finite e the change in F has a very complicated 

dependence on e. But, since e is assumed small, we can Taylor expand with 
respect to and then "throw out" all terms beyond the term linear in eO. 
This makes sense since we assume that e is a small as we want. Thus, the 
linear part of 8F is a functional of and it is a linear functional. We denote 

SF 

the linear response kernel of this functional by the symbol Sn( -y called "the 
functional derivative". Therefore we can write: 

f SF 

SF = F[n + eO] - F[n] = ——e(S>(f)d 3 r + 0(e 2 ) (1.10.7) 

J Sn(r) 

SF 

Note, that this the linear response kernel is independent of O(r) and e. 

One technique to obtain the kernel is taking O(r) = S ? (r) = S(r — F). To make 
this notion precise, we define: 

SF F[n + eS f ]- F[n] „« nm 
- lim— — (1.10.8) 



Sn(r) e^o e 

Where S ? = S(r — F) is a delta function, expressing the fact that we make a 
localized change in n(r) at F. 
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Examples 

Apply definition (1.10.8) to F rQ [n] , and I v [n]. Find the functional derivatives 
and check that Eq. (1.10.7) holds. 

Solution: 

8 F ro ,. F r<\ n + eS f] ~ F r [n] (n(r ) + eS f (r )) - n(r ) 

= hm = km- 



£n(r) £^o e e^o e (1.10.9) 

= 8 f (r Q ) = S(r Q -r) 

The functional derivative comes out S f (r) = <5(r — f), meaning that any 
change in n(r) made at r will not affect /v„[ n ]/ unless it is made at f = r . 
Indeed, if we make an arbitrary small change n -> eO, we have from (1.10.7): 

SF ro = F rQ [n + eS f ] - F rQ [n] = (n(r„) + eS f (r )) - n(r ) = 5 f (r ) 
= 5(r - f) 

And this is seen to be indeed the case directly from the definition of F r in Eq. 
(1.10.1). 

As for I v [n] =-f v n(r)d 3 r : 



(1.10.10) 



Sl v I v [n + eS f ] - I v [n] 
~ hm 



Sn(r) e^o 



= fy«r) + eS- r (r))d 3 r - V' 1 f v n(r)d 3 r (1 10 n) 

V' 1 j 8 f (r) d 3 r = V' 1 



Another Example 

An important functional, we will use is the "Hartree energy", the classical 
repulsion energy associated with a charge distribution n(r): 

Mn]=i//^^V (1.10.12) 
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Prove that the functional derivative of the Hartree Energy, called the Hartree 
potential is the electrostatic potential of that density of electrons: 

SE H 



v H (r) 



Sn(r) 



1 / (n(r) + eS f (r))(n(r') + eS- r (rO) _ n (r)n(r') \ 3 3 
2 JJ \ \r-r'\ \r-r'\ 



= lim- 

e->0 



1 (eSr{r)n{r')+n{r)eS f (r/)) , 

; — - '-d^rd^r 

.. 2 JJ r-r' 

= km 



(1.10.13) 



e->0 

J |r-r'| 



Thus, we find as requested that v H (r) is ineed the electrostatic potential of the 
density n(r). 

Two shortcuts: when "functional integrating", one can often use regular 
differentiation rules. All one needs to remember is the following shortcuts: 

———^8(x-x') ^>8{x-x') (1.10.14) 

8n{x') 8n{x') 

And use them within "chain rules": S #^P ~» f'(n(x)^)8(x — x') (where 

f'(n) = ^~/( n ))- An example of the use of this rule is the previous exercise: 

SI C 8yi(y) C 

= V- 1 ^r-r4d 3 r = V- 1 8(r - f) d 3 r = V' 1 (1.10.15) 
8n(r) J 8n(r) J 

Another example is the following functional, related to the kinetic energy: 

p CO 

Jlf] = f'ixYdx (1.10.16) 

J — 00 

The functional derivative is: 
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(1.10.17) 



= J 2f'(x)8'(pc - x)dx = -2f"(x) 

J — CO 



Here we used the identity: f g(x)8'(x — x)dx = — g'(x) and the shortcut 
= S'(x - x) 

8 fix) OK - X X) - 

iv. Invertible function functionals and their functional derivatives 

Consider a "continuous set of functionals" which is "indexed" by a 
continuous variable r in some domain D. More precisely u[n](r) assigns a 
functional to each r in fl. Let us suppose further that n is a function of r 
defined on the same domain D. An example is the relation between the 
electrostatic potential and the charge density: 

v H [n](r) = j j^\d 3 r' (1.10.18) 

For any charge distribution n(r') the integral forms a the electrostatic 
potential v H (r). The potential at a give point r depends in principle on the 
charge density at all points in space. This is a very non-local dependence. 

What about the functional derivative? If one makes a small change in the 
function n at the point r', 8n(r'), how will that affect the functional at r? We 
write this as a definition of the functional derivative: 

8v[n](r) = —LpJ-Sn(r')d 3 r' (1.10.19) 
J 8n(r') 



In the example above, it is immediately noticeable that: 

8v H [n](r) 1 
8n(r') \r — r'\ 



(1.10.20) 
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An interesting concept can arise when one knows in advance that the relation 
v[n] is unique, i.e. 

{v[n](r) = v[n'](r) for all r} <=> {n(r) = n' (r) for all r} (1.10.21) 
In this case one can also consider n to be a functional of v. Indeed: 

Sn[v](r) = / ;\ J Sv(r')d 3 r' (1.10.22) 
J ov(r ) 

Once can combine this equation with (1.10.20) and see that: 

fcM = f *WW f „ v (1 . 10 .23, 

v J J Sv(r') J Sn(r") v } K ' 



This shows that: 



/ 



Sn(r) Sv(r') „ 

— HV- ; , d 3 r f = S(r - r") (1.10.24) 
5v(r')5n(r") 

Thus, there is an inverse relation between the functional derivatives. This is 
similar to what happens in usual derivatives with inverse functions. When a 
function g{x) is invertible we can speak of x(g). The change in g due to a 
change in x is: dg(x)=—dx. Similarly: dx = —dg. Therefore: ( — ) = 



dg 

, -1> 



v. Convex functionals 

A function f{x) is said to be convex if for any pair of abscisas x x and x 2 the 
curve (x,/(x)) described by f(x) in the interval x 1 <x<x 2 is always beZow 
the straight line connecting (xi,/(xi)) and (x 2 ,/(x 2 )). A more formal 
definition is that for any < A < 1 we have: 

f(Ax 2 + (1 - AK) < Af(x 2 ) + (1 - A)/(Xi) (1.10.25) 
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A useful interpretation of the above definition for convexity is in terms of 
averages. We can view Xx 1 + (1 — X)x 2 as a weighted average of the two 
points, with positive weights (such an average is usually called a convex 
average). In this sense we have convexity obey: 

/«*)) < {fix)) (1.10.26) 

This relation is much more general than just 2 points . We can easily show that 
for convex functions it holds for any averaging procedure (x) — £t QXj and 
(TOO) = Y,i c if(. x i) with c £ > and £i c i = 1- Eq. (1.10.26) is sometimes called 
Jensen's inequality. 

One of the important implications of a convex function is that it cannot have a 
local minimum: either there is no minimum (for example f{x) = e x ) or there 
is just one global minimum. If a function is known to be convex and to have a 
minimum then any "descent" algorithm is surely to find the (global) 
minimum. This is useful. 

We can Taylor expand f{Xx 2 + (1 — X)x 1 ) = f{X{x 2 — x x ) + x x ) with respect to 
X, assuming now X « 1. Then: 

Xf(x 2 ) + (1 - X)f{x 1 ) > f{X{x 2 - x x ) + x x ) 

(1.10.27) 

= + V/Oi) • A(x 2 - *0 + 0{X 2 ) 

After rearrangement and division by X we set X to zero and obtain: 

/(x 2 ) - f{x t ) > V/Cq) • (x 2 - Xl ) (1.10.28) 

This relation, if obeyed for all pairs of points, is equivalent to the convexity of 
the function. By exchanging x x and x 2 we obtain also: 

V/(x 2 ) • (x 2 - Xl ) > /(x 2 ) - f(x x ) > V/(Xi) • (x 2 - x x ) (1.10.29) 

From (1.10.28) we further have, for A -> 0: 
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f(x x ) + IV/ (x x ) • u < /(Xi + lit) 

1 

= /(xj + AV/CxJ • 8x 1 + A 2 - 8x ± ■ VV/CxJ • 8x 1 (1.10.30) 
+ 0(A 3 |u| 3 ) 

From this, after cancelling and dividing by X we find -u ■ VV/Xx-l) • u > for 
arbitrary u which means the Hessian VV/Xx-l) is positive definite: another 
characteristic of convex functions 

A similar definition of convexity can be applied to functionals 

F[An x + (1 - A)n 2 ] < AF[nJ + (1 - l)F[n 2 ] (1.10.31) 

And they too have the equivalent differential property which will be useful 
below: 



F[nJ - F[n 2 ] > j (^^y) (%(r) - n 2 (r))d 3 r (1.10.32) 



5F 

An example for a convex functional useful for density functional theory is the 
so-called Hartree energy functional: 

1 ff n(r)n(r') 



■d 3 rd 3 r' (1.10.33) 



\r — r'\ 

The second functional derivative of this function is |r — r'| _1 which is a 
manifestly positive matrix (this is easily seen by Fourier transforming the 
function v(r) = with result v(q) = which is positive). Another example is 

the von-Weizsacker kinetic energy functional: 

E vW [n] = -i| ^Mr)V 2 ^Mr)d 3 r (1.10.34) 
The physical meaning of this functional is that it gives N times the kinetic 

jn(r) 

energy of a single-particle with ground state wave function ip(r) = 1^— 
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where N = / n(r)d 3 r. Again, by showing the Hessian is positive definite we 
can straightforwardly show this functional is convex First, rewrite it as: 

E vW [n] = i| (vV<0) 2 d*r = ±j ^d 3 r (1-10.35) 

We drop the explicit r in the last expression, for brevity Then, make an 
arbitray but small perturbation eO(r): 



1 1 f (Vn + 6rVO) z _ 

[n + eO] = - d s r 

8 J n + 



(1.10.36) 



And expand the denominator to second order in e: (n + eO) 1 -> ^1 — e* + 

e 2 + "-^ . Expanding the numerator and performing the multiplication 

order by order one can then read off the second order contribution and after 
some manipulation bring it to the form: 



e 



^E vw [n + e^>] = € -j -(^" V) d3r (L1037) 



Clearly, this second order term is absolutely non-negative for all perturbations 
and at all physical densities. This shows that the underlying functional 
derivative Hessian is positive definite and thus the functional is convex. 

A similar condition exists for the second functional derivative. The Jensen 
inequality holds here as well: 

F[(n)] < (F[n]) (1.10.38) 

This relation, in combination with the fact that f{x) = e x is convex is useful 
for developing mean field approaches in Statistical mechanics. It was used by 
Gibbs and later, in a quantum context by Peierls.[9] 



Electron Density Functional Theory 
© Roi Baer 



Page 57 



J. Minimizing functionals is an important part of 
developing a DFT 

Equation Section (Next)Minimization of functions 

We start with the problem of minimization of a simple ID function. Given a 
function fix), we want to find a point x* such that the function is minimal. It 
is clear that the slope of fix) at x = x* must be zero. If the slope is positive 
then we can go left (decrease x*) and reduce fix). Same logic (but to right) if 
it was negative. Thus, a necessary condition for a minimum is: 

f'(x*) = (1.11.1) 

Let us expand in a Taylor's series the function around the point x*. Clearly we 
have: 

fix) = f{x*) + f'ix*)ix - x*) + i/"0*)0 - x*) 2 + ... (1.11.2) 

However, taking Eq. (1.11.1) into account, 

fix) = fix*) + - x*) 2 + - (l.H.3) 

When x is extremely close to x* we may neglect high order terms and then we 
see that Eq. (1.11.3) is the equation of a parabola. In order that fix) be a 
minimum, we must have an ascent of fix) when we move away from x* , no 
matter which direction we take (left or right). For small displacements x* -> x, 
we see from Eq. (1.11.3) that the change in fix) is -f"ix)ix — x*) 2 . Since 
(x — x*) 2 is always positive, no matter the direction we move, we must 
demand that fix*) > as well. Thus, our necessary conditions for a 
minimum are: 

/'(**) = and fix*) > (1.11.4) 
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Now, let us consider the case of functions of two variables, f(r) where 

(X\ 
yj = (x y) r . Notice that we use the "transpose" symbol T to turn a 

column vector into a row vector. Let us Taylor expand f(r) around a point r*\ 
fir) = f(r*) + V/(r*) T (r - r*) + -(r - r*) T VV/(r*)(r - r*) + - 



(1.11.5) 



= f(r*) + g* T (r -r*) + -(r- r*) T H*(r - r*) + 



Note that we use the notations: 



g* = g(r*) = V/(r*) = 




H* = H(r*) = 



1 

Sx 2 SySx 
S 2 f S 2 f 



(1.11.6) 



\8x8y Sy 2 ] , 



The vector g{r) is called the gradient and the matrix H(r) the Hessian at r. 
Note that the Hessian is a symmetric matrix, since the order of differentiation 
is immaterial. When r* is a minimum, moving infinitesimally away from it in 
any direction will not change the function. Why? We can show this leads to a 
contradiction. Let d be an arbitrary direction. If the function decreases when 
moving from r* in that direction, no matter how small the step size, then this 
contradicts that r* is a minimum. Now, if, d is an ascent direction, then — d 
will be a descent direction - again a contradiction to r* being a minimum. 
Thus we conclude , that / cannot change (to first order) or in other words, the 
gradient g* of f(r) at r* must be zero. Then Eq. (1.11.5) becomes: 



f(r) = fir*) + -(r - r*) T H\r - r*) + 



(1.11.7) 
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Continuing, in order for this function to have a minimum at r* the second 
term on the right hand side must always be positive. Thus Hessian matrix 
must be such that for any vector v T Hv > 0. A matrix having this property 
is called "positive definite" (PD). When a PD matrix is symmetric, then its 
eigenvalues must be strictly positive. We can summarize the necessary 
conditions for a minimum at r*: 



As an example, let us take the function f{x) = x 2 + y 2 . This function is always 
non-negative, so its minimum is at x = y = 0, where it obtains the value 0- 
The gradient is indeed zero at (x, y) = (0,0): 



And the Hessian is 2I 2 where I 2 is the 2><2 identity matrix, thus it is trivially 



ii. Constrained minimization: Lagrange multipliers 

Next we discuss constrained minimization. Suppose we want to minimize 
f(x>y)/ under a constraint that y = g(x). This is "easy", since it can be 
transformed into an unconstrained minimization. All we have to do is 



flf * = V/(r*) = and WisPos.Def. 



(1.11.8) 




(1.11.9) 




PD. 



minimize F(x) = f(x, g(x)). This will give: 




(1.11.10) 
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+ 



'd 2 f S 



+ ^ (dxdy) 



g'ix) + C£) g"{x) 



(x,g{x)) 

And we can then search for x* such that F'(x*) = and F"(x*) > 0. 

But sometimes it is not possible to write the constraint directly as y = g(x). A 
more general form is: h{x, y) = 0. This is also more symmetric. Consider the 
situation of minimizing f(x,y) under the constraint /i(x,y) = 0. Suppose the 
point depicted in the following figure is this minimum. 




h=l 



h=0 



h=-l 



This means that moving slightly along the isopleth of h = from r* will not 
change /: if it decreases / this is not a minimum and if it increases / we will 
move in the opposite direction and / will decrease. 

Thus at r* the gradient of f , V/(r*) 7 must be normal to the line h = 0. 
Furthermore, by definition the gradient of h is normal to its isopleth. Thus, the 
necessary condition for minimum is that both vectors point to the same 
direction and so are proportional. We denote the proportionality constant by 

A*: 
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V/(r*) = A*Vh(r*) 

(1.11.11) 

h(r*) = 

Lagrange noticed that this relation and wrote a Lagrangian for the constrained 
problem: 

L(r,X) = f(r) - Aft(r) (1.11.12) 

The new variable 1 is called a "Lagrange multiplier". The necessary conditions 
for a minimum in Eq. (1.11.11) translate into the condition of finding a 
stationary point of L 

VL(r*,A*) = (1.11.13) 

Where V= (d x , d y , d z , d x ). The derivative with respect to A gives us back the 
constraint. As an example, consider minimizing /(x,y) = x + y under the 
constraint /i(x,y) = x 2 + y 2 — 1 (i.e. find the point on the unit circle for which 
x + y is minimal). We have: 

L(r, X) = (x + y) - A(x 2 + y 2 - 1) (1.11.14) 

The necessary condition is: 

VL = (l - 2/Tx*, 1 - 2X*y*, -(x* 2 + y* 2 - l)) = (0,0,0) 



1 / 1 * 2 



(1.11.15) 



X _y ~2l*' V2A* 



The solution is thus: 



1 1 

A* = + — x *=y* = +— (1.11.16) 

V2 V2 

Since these are necessary conditions, we need to consider them further (they 
might correspond to a maximum). The minimum corresponds only to the 
negative solution. 
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The method of Lagrange multipliers is often more convenient to work with 
than direct replacement. The problem is thus that of finding a minimum of 
f(Xi, under the M < N constraints ft(r 1( - ,r N ) = c, or /ij(r 1; ...,r N ) = c t 

i = 1, ...,M is: 

find r* such that: h(r*) = c for which : f(r*) < f(r) 

(1.11.17) 

for any r such that: h(r) = c 
To facilitate such a search, we formulate the Lagrangian function: 

M 

L{r,k;c)=f{r)-Y j h(h i (r)-cd (1.11.18) 

t=i 

Our plan is to find the position of r which minimes L for any choice of A and 
then change c until = c. At A* we have L assuming a minimum at a point 
r* . A necessary condition for the constrained minimum to be achieved at the 
point r* and with the Lagrange multipliers A* is: 

V r L(r*,r,-c) = 

(1.11.19) 

V x L(r*,r-,c) = h(r*)-c = 

Note that A* is not necessarily a minimizer of L. In fact the opposite is tue: A* 
is a maximizer of L and the search for constrained minimizations is a search 
for saddle points (see ref. [10] for a method to solve such a problem on the 
large scale). 

It is interesting now to ask how f(r) changes if we change the value of the 
constraint Cj. Indeed, when the constraints are changed, the optimized point 
and Lagrange multiplier can change, so the Lagrangian is changed: 

5L(r*,r; c) = L(r* + Sr*,A* + SA; c + 8c) - L{r*,X*; c) (1.11.20) 

Now, because of Eq. (1.11.19) we find: 
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SL(r*,X*; c) = L(r* + Sr*,X* + SX*; c + 8c) - L(r*,X*; c) 

(1.11.21) 

= L(r*,X*; c + Sc)- L(r*,X*; c) 



Thus: 



d , ^5 
—L(r*,X*;c) = - 
aci oCj 



—L(r*,X*;c) = —L(r*,X*;c) = X\ (1.11.22) 



Since L(r*,X*; c) = f(r*; c) we find: 

^-/(r*;c)=AJ or V c /(r*; c) = A* (1.11.23) 

This equation reveals the "meaning" of the Lagrange multipliers X* at the 
optimal point: they are equal to the rate at which the optimal value of the 
minimized function / changes when c ir the value of the i constraint, is 
changed. This is an important result which we use below whenever we want 
to give physical significance to Lagrange multipliers. 

iii. Minimization of functionals 

The same considerations for functions apply for functionals. Given a 
functional / (/), a necessary condition for its minimum is: 

SI 

= (1.11.24) 



8f(r) 



For example, consider a ID classical particle of mass 77 in a potential well 
v(x). The action S is a functional of the trajectory of this particle: 



S[x] 



-r 



— mx(t) 2 + v{x(t)) 



dt 



(1.11.25) 



For any trajectory x(t) between times t and tp L[x] returns a number. 
Lagrange showed that finding the trajectory that makes L[x] stationary 
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(although, not necessarily minimal), under the condition that x(t ) = x and 
x(tf) = Xf are given and thus not varied, is equivalent to solving Newton's 
equations under these same constraints. The functional differentiation of the 
kinetic energy is performed with: 



- x(t') 2 dt' = 2 x(t')-vrdt' = 2 x(t')5(t - t')df 
= -2 J x(t')5(t - t')dt' = -2x(t) 



The condition for stationarity of the action under changes in the trajectory 
which leave the edges intact give: 

SS 

= = -mx(t) - i/(x(t)) (1.11.27) 

From which we obtain Newton's equation mx(t) = — i/(x(t)J. This equation 
must be solved subject to the given constraints at the endpoints, i.e. x(t ) = x 
and x(tf) = Xf. 

Minimizing functionals with constraints can again be done using Lagrange 
multipliers. Then one defines: 

L[f,X\=I[f]+M[f] (1.11.28) 
And the necessary conditions are: 



SL 



Sf(r) 



SL 

= M 



= h[f*] = (1.11.29) 



For example, let us solve the following problem. Find the shape in a plane of a 
closed contour encircling maximal area under the constraint of a given 
circumference l . We need to find the contour given by r(6>) with r{B + 2n) = 
r(0) which gives zero variation to 
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L [r, X] = S[r] + A(L[r] - Z ) (1.11.30) 

Where S is the area and L is the circumference. We limit ourselves to curves 
that can be uniquely defined. A point on the curve, is defined by the angle 
and the distance from the origin r(0). The vector to a point on the curve and 
its derivatibe is: 

r(0) = r(0)(cos0,sin0) 

(1.11.31) 

r'(0) = r'(0)(cos 9 , sin 0) + r(0)(-sin , cos 0). 
The area of a small arc from to + d0 is dS = - |r(0) x r(0 + d6)\ = 
-|r(0) x r'(0)|d0. Notice that |r x r'| = r 2 and so d5 = ^r(0) 2 d0. 

2 2 

The square circumference is: dL 2 = (r(0 + d9) — r(0)) = (r'(0)) d0 2 . Thus 
dL = \r'{9)\d9 = Vr(0) 2 +r'(0) 2 d0. Thus, the area is a functional is: 

S[r] = J -r(0) 2 d0 

(1.11.32) 

/-27T 

L[r] = V r (0) 2 +r'(0) 2 d0 
•'o 

The functional derivative of S is easy: 

S^-rW (1-11.33) 
The functional derivative of L is computed by adding a "small" function e(0): 



<5/ = J'' ^(r + e) + [r' + e'JdO' - 1 



= (* 4r 2 + 2re + r n + 2r'e'd6' - I (1.11.34) 

*J 



2^ + 2^ 

2 . 12 

r + r 



Then, we can use yji + n = i + -r? + (r/ 2 ) so: 
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61= f 2T r f + r£ d9' (1.11.35) 

JO ^2 , _,2 



r 2 + r /2 



Integrating by parts the second term (remembering that the end terms drop 
because G — 2n is the same point as 6 = 0: 



r z " re r" / r y 

= — d0' - — ) edO' (1.11.36) 

J Vr 2 +r' 2 J Wr 2 +r' 2 / 

We thus get: 



5Z 



(1.11.37) 



<5r(0) ^V 2 _|_ r '2 Wr 2 + r' 2 ' 
We define the Lagrangian: 

L[r]=5[r]-A(I[r]-I ) (1.11.38) 
And get the maximum from: 

SL ( r ( r \'\ 

= — — = r - I ( - ■ (1.11.39) 

ffr(0) Wr 2 + r' 2 Wr 2 +r' 2 / / 

Clearly, a solution is r(6>) = 1, a constant, since then all derivatives are zero. 
The curve is then circle of radius X and given the circumference, we have: 
X = l /2n. 

Exercises for chapter 1 

1) Consider a system of two particles in a harmonic well. Each particle repels the other. 
The potential energy is V (x 1 , x 2 ) = ^ kx\ + ^ kx\ — ^ q — x 2 ) 2 . We assume that 
the system is bound (so < 2q < k). (a) Write down the exact Hamiltonian for this 
system, assuming the mass are m 1 and m 2 . (b) Make a transformation separating 
the Hamiltonian into two uncoupled particles: the center of mass X = 1711X1+7112X2 

m 1 +m 2 

and the reduced mass x = x 2 — x x . Write the Hamiltonian of each part and 
determine the mass of each particle. What are the masses when m 1 » m 2 l (c) 
Write down the energy levels of the system and the eigenfunctions. (d) Now make a 
Born-Oppenheimer approximation. Assuming m 1 » m 2 , determine the Born 
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Oppenheimer potential surface (obtained when x x is considered a parameter). 
Solution: The Hamiltonian is: 

ft 2 , ^ , 1 , 1 , 1 , 
= -- — d} -- — di + -kx? +-kxl --q{x 1 -x 2 ) 2 . 
2m 1 1 2mi 2 2 1 2 2 2 WV 1 2y 

In the new coordinates: 

(m 1 x 1 + m 2 x 2 ) m 2 m 1 

X — : , X — X-i X-\ , X-i — X z~i~X. Xn — X H X. 

M M M 

If /(x 1( x 2 ) is any function than we can define the "same" scalar function as 

F(.x,X)=f{x-Jx,X + ^x). 
2) Calculate the correlation in the case of a singlet, where both electrons are in the 
ground state orbital and when one is in the ground state orbital and the other in the 
first excited orbital. 



II . My first density functional : 
Thomas -Fermi Theory 

A.Basic concepts in the electron gas and the 
Thomas -Fermi Theory 

Thomas and Fermi assumed that the energy of an atom or a molecule can be 
written as a functional of the 1-particle density as follows: 

E TP [n) = T TF [n] + f v ext {r)n{r)d*r + d 3 rd 3 r' (2.1.1) 

(Note, for use via the Born-Oppenheimer approximation, to this energy we 
need to add the nuclear-nuclear repulsion energy) They then assumed that 
the density that characterizes the ground-state minimizes this functional 
under the constraint: 

j n(r)d 3 r = N e (2.1.2) 

The first question, beyond the rigor of this approach is, what is the kinetic 
energy functional? In order to take into account the Fermi nature and the 
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quantum nature of the electrons, this functional must include both these 
considerations. The Thomas Fermi solution is to assume: 



What shall we take for t s (ri)? Consider first a simple case: a homogeneous gas 
of density n (i.e. n(r) is independent of r). Furthermore, let us assume that the 
electrons are non-interacting. This is a simple enough system to enable the 
analytic calculation of the kinetic energy functional. From the form of (2.1.3) 
we see that the total kinetic energy is the sum of contributions of various 
infinitesimal cells in space. Each cell contains n(r)d 3 r electrons and so, if we 
interpret t(n) as the kinetic energy per electron of a homogeneous gas of non- 
interacting electrons then this sum is yields exactly the total kinetic energy for 
this homogeneous gas. The Thomas-Fermi approximation then uses this same 
t(n) also for the inhomogeneous interacting case. 

Let us now compute t(n). Consider a homogeneous gas ofN uncharged electrons. 
They are non-interacting. These electrons are put in a cubic cell of length L. The 
electron density is everywhere the same n = — = 



We assume the wave functions are periodic in the box. According to Fourier's 
theorem, we can write any periodic wave function as a linear combination of 
plane-waves, as follows: 



T[n] = f t s (n(r))n(r)d 3 r 



(2.1.3) 



v 




(2.1.4) 



n 



Where: 



1 [Jx> ly> ^z) 1 



(2.1.5) 
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and l x /y/z are integers. Fourier's theorem is based on the orthonormality of 
the plane waves 

GftdlM = \ fff e £ ( k - k '> r d 3 r = 5 k y (2.1.6) 
Where we defined 

e ikr 2n 

Wr) =VF k = T" <2 ' 17) 

We imagine 3-dimensional k-space divided into an array of small 
compartments, indexed by a set of integers n = (l x , l y , Z z ) or by the vector k. 

"2.71 ( 2tz0*^ 

Each compartment is of k-length A/c = — and its k- volume is A/c 3 = - . For 
large r-space boxes the k-space compartment is extremely small since A/c 3 is 
proportional to the inverse box volume. Since we are interested eventually in 
the limit L -> oo, we may assume approximate sums of any function /(k) over 
the discrete values of k = — n by integrals: 

^/(k)^^J7(k)d 3 /c (V^oo) ( 2.i.8) 

k 

Let's show that plane- waves are eigenstates of kinetic energy operator 7^: 

h 2 1 h 2 k 2 

Now, consider the wavefunction of the N e non-interacting electrons in their 
ground-state. Since they are non-interacting, this wave-function is a product 
of single-electron wave-functions: 

V = iMrO^fod -^k Ne/2 (r We - 1 )^ kNe/2 (r We ) (2.1.10) 

Here i/^k( r ) is the state of a spin-up electron with wave vector k. while i/^k( r ) is 
the state of a spin down electron with wave vector k. Anticipating the 
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antisymmetry, we build this wave function by placing 2 electrons in the same 
spatial orbital (once with spin up and the other with spin down). Since non- 
interacting electrons have only one type of energy, i.e. kinetic energy: 

H = Zn=i (~T~^n) r we can easily show that (2.1.10) is an eigenstate of the 

Hamiltonian: 

Z N e ( h 2 \ 
n __\-^) K to) ' * ' ^k We/2 (r We -i)^ kNe/2 (r w J 

/h 2 k 2 \ 

= 2, n=1 2 (^fj^Cri)^^) ■■'^ Ne/2 ^N e -i)^ Ne/2 ^N e ) (2-1.11) 

n=l V 2m e/ 

One sees that the energy is just the sum of kinetic energy E n =i 2 ( — -) in each 

spin-orbital of the product wave function. Let us now antisymmetrized this 
product wave function. We do this by adding all products resulting from even 
permutations of the electrons and subtracting all odd permutation products. 
One convenient way to represent such a sum is using a determinant, called a 
Slater wave function: 



¥ = 



/JV! 



:det 



fa) 

0k! ( r H<L 



(2.1.12) 



For this wave function to be minimal energy must fill 2 electrons per level 
starting from the lowest kinetic energy and going up until electrons are 
exhausted. Denote the highest filled level by k F . Then: 

V 



N, 



filled 



(2n) 



■jjj 9(k F -k)d 3 k. 



(2.1.13) 
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Where 9{x) is if x is negative and 1 otherwise. This is called the Heaviside 
function. We now perform the integral using spherical coordinates: 

V f kp V An V kl 

= (2nr} 4Mk = (2^T k ' = SPT <2 - U4) 

The number of filled orbitals is the product of the real-sapce volum V and the 
k-space occupied state volume, divided by (2n) 3 . Since N e = 2Nf iUed and the 
density is n = — we have: 



n = 2N f illed = (2.1.15) 



V 3n 2 ' 

The electron density determined directly the highest filled momentum state: 
2V f kF fh 2 k 2 \ , 2V h 2 kl 2V h 2 kl 

™ = wl {^r kdk = w An ^i = ^^i- (2 - 116) 

The energy per particle is: 



r s l h 2 kl 3h 2 k 2 3 h 2 2 

t s (n) = T7 = -^-t. r = -p-z = — (3?r 2 n)3 = C TF n 2 ' 3 , (2.1.17) 

N n^nzrrig 5 5 2m e 5 2m e 



where: 



3 h 2 23 2 

C TF = -- — (3Ti 2 )i = — (3Ti 2 )iau = 2.871 au . (2.1.18) 
5 2m e 10 

Plugging into Eq. (2.1.3), the Thomas-Fermi kinetic energy functional is 
obtained to be used in Eq. (2.1.1): 

T s [n] = C TF f n(r)id 3 r. (2.1.19) 



Exercise: The Thomas Fermi functional for the hydrogen atom. 

a. Repeat the calculation above but now for a "spin-polarized HEG". That is, do not 
assume that there are 2 electrons in each k-state (the "spin-unpolarized" case) but 
instead, that all spins are up and so there is only one electron per k-state. 
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b. Since the electron in a hydrogen-like atom is spin-polarized, use the Thomas-Fermi 
KE functional derived in (a) and compare its estimation of the kinetic energy of the 
electron in a hydrogen-like atom to the exact value. Using the exact kinetic energy in 
the hydrogen atom (you can find it using the virial theorem), assess the quality of 
the result as a function of the nucleus charge Z. 



B. Minimization of the Thomas -Fermi energy 

Now, according to the TF theory, the true ground-state electron density is the 
one that minimizes E TF [n]. But the electron density must also account for the 
required number of electrons of N e , so there is a constraint for the 
minimization: 



j n(r)d 3 r = N e 
Thus, we must build a Lagrangian to be minimized as: 



(2.1.20) 



L[n,X\ = E TF [n] - fi 



I" 



(r)d 3 r - N e 



(2.1.21) 



Minimizing it gives the Thomas-Fermi equation: 

SL 8E TF 

= SMF) = 5MF)- fl (ZL22) 

We see that the Lagrange constant [i is the chemical potential, since it is equal 
to the change in energy when we perturb the density and this change is 
everywhere constant. We must now compute the functional derivatives: 

8 f „ .„ . 5 



8n(r) 



j C TF n(r')V 3 d 3 r' = -C TF n(r) 2 / 3 
v(r')n(r')d 3 r'^j = v(r) 



8n(r) 
8 1 ffn(r")n(r') 



(2.1.23) 



n(r') 



d 3 r' 



8n(r)2JJ \r" — r'\ J \r — r'\ 

Plugging into Eq.(2.1.22), one obtains the Thomas-Fermi equation for an atom: 
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(i = |c rF n(r) 2 / 3 + v(r) + [ , n(r) „ dV (2.1.24) 
3 J \r — r'\ 

This is an integral equation for n(r). It is called the integral Thomas-Fermi 
equation. The potential v(r) is due to the positive charge, hence we can write: 

n fr'l 

V (r) = -J^-ZdV.Now we can define a fofaZ potential 



, n.(r') — n(r') , 



as the sum of the total electrostatic potential and the chemical potential. Since 
V 2 - = AnS{r), we have: 

r 

-V 2 0(r) = 47r(n+(r) - n(r)), (2.1.26) 
leading to an equation for the total potential: 



1 _ 

— V 2 = 

47T 



3 



0(r) 



3/2 



n+ (r) (2.1.27) 



.5C TF 

This is called the "differential Thomas Fermi equation". The constant fi is 
buried in but it did not disappear: it must be chosen so that Eq. (2.1.20) is 
obeyed. Furthermore, it is clear that the potential 0(r) is manifestly non- 
negative in TF theory (otherwise we could not take its square root). 

As we showed (see Eq. (1.11.23)), [x can be shown to act as a chemical 
potential, i.e. -— = mOV) We will not solve this equation for atoms or 

molecules. We just comment that it gives a smoothed value for the atomic 
density, not showing the shell structure. We can see this in the following 
figure, where we plot the radial density as computed by a relatively accurate 
theory, such as Hartree-Fock and the TFD (Thomas-Fermi-Dirac) theory. 
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Kinetic Energy poorly 
represented! 



1.0 



Sqrt[r] a.u. 



2.0 



There is a question of how does the minimal energy of the Thomas Fermi 
functional compare with the accurate quantum mechanical energy. This 
question has been examined. It was found that for atoms with Z -» oowe have: 



E TF {N = XZ) 
lim ; T=r = 1 



(2.1.28) 



z-*™E exact (N = XZ) 

For < X < 1 (i.e. the number of electrons is smaller than that of the protons 



N 



and - is held while Z -» oo). Note that the Thomas Fermi energy for an atom 



has the property that: 



Eatomi^-Z) — Z 7 / 3 Ea[ om (_X,l) 



(2.1.29) 



C. Thomas-Fermi does not account for molecules 

Consider an existing system, and now add a bit of positive charge 8q = 
j 8n + (r)d 3 r, where 8n + (r) > everywhere. Also, increase the electronic charge 
8q so as to preserve neutrality. The electron distribution is changed by 8n(r) 
(so that 8q = J 8n(r)d 3 r). The change in the total energy 5£"[n] is composed 
of the change in electronic energy fiSq, and the corresponding change in the 
electrostatic positive charge energy: 
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8E[n] = 8E TF [n] + 8E NN [n] 



(n + (r) — n(r))Sn + (r') 
\r — r'\ 



d 3 rd 3 r' 



(2.1.30) 



liSq + 



Using the definition of the total potential in Eq. (2.1.25) we find: 



5E[n] = I (p(r)Sn + (r)d 3 r 



(2.1.31) 



Since everywhere we add positive charge, i.e. 8n + (r) > and since 0(r) is non- 
negative everywhere, this process increases the total energy! In TF theory 
addition of infinitesimal positive charge followed by addition of a 
compensating electronic charge causes an increase in total energy of the 
system. In a more elaborate treatment Teller showed[ll] that the total energy 
of any diatomic molecule is higher than the sum of energy of its constituent 
atoms, i.e. that TF theory cannot account for stable molecules! He concluded 
that Thomas Fermi theory is not very useful for chemistry. 

D. Thomas-Fermi Screening 

When a point impurity Ze is inserted into an electronic system, it pulls (Z 
positive) or repels (Z negative) electrons towards it. This has an effect that the 
impurity is partially screened by opposite charge and so it has a smaller effect 
on distant charges. Let us study this phenomenon in the electron gas, using 
Thomas-Fermi theory. The homogeneous gas of electrons is a model for ideal 
metals, so the screening effect we address here is relevant for many metallic 
systems. Macroscopically, the "free" metal electrons completely screen the 
charged impurity. However microscopically, perfect screening is not possible 
because electrons have kinetic energy - even at zero temperature - and a short 
ranged electric field develops around the impurity. Thomas Fermi theory 
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takes kinetic energy effects into account and can be used to estimate the form 
of the local electric field, specifically its size or length scale. 

Let us study an unperturbed homogeneous electron gas using Thomas-Fermi 
theory. Such a "gas" has no structure and it is characterized by only one 
parameter: its density n . In order to neutralize it and support the electron 
homogeneity we add positive smeared homogeneous charge density +en . 
All the Coulomb energies (e-e, e-N and N-N) cancel exactly so the only energy 
left is the electronic kinetic energy: 

E TF [n ] = j C rK n (r) 5 / 3 d 3 r (2.2.1) 

The constraint minimization of this functional yields the following condition, 
relating the density to the chemical potential: 

Comparing with Eq. (2.1.15), and using Eq. (2.1.18) we find for the chemical 
potential: 

_8E n __S (k 3 F \ 2/3 _h 2 k 2 F 

Thus we see that indeed the electron density is constant and the chemical 
potential is equal to the kinetic energy corresponding to the maximal 
occupied momentum hk F . 

Now we introduce a positive charge Ze. The density of electrons is changed: 

n(r) = n + n^r) (2.2.4) 

It is physically clear that n 1 (r) is localized around the impurity (assumed at 
the origin). We therefore have for the total energy of the system in terms of %: 
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(2.2.5) 



2 JJ \r-r'\ 

The corresponding TF equation comes from minimizing: 

M = ^ = ^ rF (n 0+ n l( r)) 2/3 -^ + ^ f^^dV (2.2.6) 
^ 57ii(r) 3 " v 1 y r J |r-r'| 

n fr') 

We write: 0(r) = e f , 1 ,, rf 3 r' and so: 

|i — r'| 

SEtp 5 , ,2/3 Ze 2 

M = 5^) = 3 ^ + ni(r) ^ T" + 60(r) (ZZ7) 

Upon linearizing, assuming n x « n : 

I c ""» /3 ( 1 + 5^)-T ! + e * (r) = '' (2 ' 2 ' 8) 

5 2/3 

We can write: - C TF n = fi and so: 

^ Ctf^q 1/3 ^i - — + e0(r) = ix - Ho (2.2.9) 
Finally since V 2 = —\nen x we have: 

-^T< 1/3 ^fV 2 - - + 0(r) = (2.2.10) 
9 47re 2 r e 

i 

We have from Eq. (2.1.15) k F = (3n 2 n )3 and we use the definition of the 

ft 2 

Bohr radius a = — — — , defining the Thomas Fermi screening parameter k TF : 

10 _-, /q 7ra n 1 

n o 7 Ctf=tt^ = t^ (2.2.11) 



9 47re 2 TF Ak F ~ /c 2 F 
With this we have the equation: 



V 2 = /c 2 F (0(r)-^-^) (2.2.12) 
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Passing to spherical coordinates we find: 

i(r0(r))" = fe| F (0(r)-^-A^). (2.2.13) 

Defining x = r we find: 

X" = kh(x-Ze-Afir) (2.2.14) 

The homogeneous equation is x" = ^tfX which has the solution Xh = 
A e -k TF r + g e k TF r Qearly, f or a localized potential solution we must take 
B = 0. To this we need add any solution of the inhomogeneous equation 
which clearly is Xm = Ze + Afir. Thus: 

X = Ae~ k TF r + Ze + Afir (2.2.15) 

This leads to: 

= — + — + Afi (2.2.16) 

r r 

In the limit that r -> we must have r0(r) -> since the electronic charge n x 
has no cusps. Thus A = —Ze . The total electrostatic potential is 

Ze —Ze 

4>tot(x) = 0(r) = e- k ^ r + An (2.2.17) 

r r 

Aside from the constant ApL, far from the impurity the surface integral of V0 tot 
evaluates to zero and by Gauss's theorem a large sphere around the impurity 
includes zero charge in it, meaning that the total amount of electronic charge 
pulled into the sphere is exactly equal to that of the impurity (Z). 

It is interesting that the screening length is proportional to k F ^ 2 or to . 

The higher the density the smaller the length, i.e. the more efficient is the 

screening, however, the dependence on n is mild because of the small 

exponent. It is also interesting to note that k TF is independent of Z. However, 

this latter results holds only in so far as our linearization is valid. For strong 
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impurities the non-linear equation will give a different result and the 
screening will depend on Z. 

E. Von Weizsacker kinetic energy 

The Thomas Fermi kinetic energy density functional is exact in the limit of 
non-interacting homogeneous gas of electrons in an infinite box. We would 
like to mention here another density functional which is exact in a certain 
limit, i.e the limit of a single electron. In this case the kinetic energy is: 

wave functions that decay to zero at r -* 

ft? 2 

oo, one can integrate by parts and obtain T = — / (Vi/;(r)) d 3 r, stressing the 

absolute positivity of kinetic energy (it cannot be zero). Finally, if xp(r) is a 
non-degenerate ground-state it can be written as ip(r) = jn(r) and so we 
obtain the kinetic energy functional of von Weizsacker: 

T vw [n] = A- J (vV<0) 2 d 3 r (2.2.18) 
Which can be written as follows, using local wave vector: 

k(r) = -—rv (2.2.19) 
2 n(r) 

So: 

h 2 k(r) 



vW 



[n] = [ W n(r)d 3 r (2.2.20) 



This functional is now used for any density, even a many electron one. The 
variation is: 



h 2 r ( (v(n(r') + eO(r'))) 2 (v(n(r')))' 



r I (v(n(r' 
J I n(r' 



ST ™ = We\ n(rO + eO(r') W)~ ' dV 
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Working this out to linear terms in e, using: (n(r') + £0(r')) 1 ~ n(r') 1 (l — 
eQi 1 nt 1 we obtain: 

Which after integration by parts of the first term finally gives: 

Thus the von-Weizsacker potential is: 

Which can be written more compactly as: 

v vW {r) = - iL [V • k(r) + k(rf ] (2.2.25) 

ZfJ. e 

Exercise: For 1-electron system, discuss the claims: 1) The wave vector fc(r) is 
the gradient of the log of the of the wavefunction: fc(r) = \7\ogxp(r) (2) the 
von Weizsacker potential is the potential for which n(r) is the ground state 
density. 

III. Many-electron wave functions 

A. The electron spin 

Zeeman has shown that a small magnetic field causes the splitting of energy 
levels in atoms. Each atomic level is split into a doublet. The amount of 
splitting is proportional to the field. At zero field these doublets are 
degenerate. The conclusion is that the electron has an intrinsic magnetic 
moment which can take two values. The states of the internal magnetic 
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moment of the electron are assumed to be proportional to an internal angular 
momentum called "spin". The spin of an electron can is assumed to have two 

values + -. This is an additional "degree of freedom". It is not continuous, but 

is nevertheless it is degree of freedom. We denote a spin-orbital ip n (x) — 
xp n (r, s) where r is a point in 3D space and s is a "spin variable", which allows 
us to perform a inner product of spin as explained now. There are two 
possible spin functions for an electron, a(s) denotes spin up and /?(s) spin 
down. These two states are complete and orthonormal: 

J = <«!«> = ! 

j /3(sTa(s)ds = (p\a) = 
j a(syp(s)ds = (a\p) = 

j P( s yp(s)ds = (p\p) = i 

The variable s is just a mneumonic. With new notation, we have: 

<^|0) = j xl)( x y<t){x)dx = j j \p(r,sy<p(r,s)dsd 3 r (3.1.1) 

B. The Paul! principle 

The electronic wavefunctions are functions of N e electronic coordinates and 
spins ^){x lt ...,x N ). Here Xj = (r ; -,5 ; ). The Pauli principle states that this 
wavefunction must be antisymmetric with respect to interchange of two 
electrons: 

i/j(...Xi, ...,Xj, ...) = -i/j(...Xj, ...,x u ...) (3.2.1) 
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This is a "boundary condition" we impose while solving for any electronic 
wave function. 

C. The Excited states of the Helium atom 

How should we represent the, in an approximate form, the low lying excited 
states of the Helium atom. He + has two low lying orbitals Is and 2s (the 2p 
orbitals are degenerate with the 2s, but we will not consider them because in 
the Helium atom they are of much higher energy. We can form a 2-electron 
wavefunctionby: xp{x 1 ,x 2 ) = ^i s (r 1 )^ ls (r 2 )[a(s 1 )/?(s 2 ) - a(s 2 )0(si)]. 

The excited states will involve excitation of an electron to the 2s orbital. We 
can then write: 

4>2(.x 1 ,x 2 ) = [xp ls (r 1 )xp 2s (r 2 ) - xp ls (.r 2 )xp 2s (r 1 )]a(s 1 )P(s 2 ) 

^ 3 (x 1 ,x 2 ) = [^isCrJ^C^) -^i s (r 2 )i/; 2s (r 1 )][a(s 1 )/?(s 2 ) 
+ a(s 2 )/?(si)] 

(3.3.1) 

^ 4 (x 1; x 2 ) = [iK s 0i>/; 2s (r 2 ) -^is(r 2 )^2 S (^i)] [/?0i)/?0 2 )] 

i/; 5 (x 1 ,x 2 ) = [iKsOi^fo) + ^i s (r 2 )i/; 2s (r 1 )][a(s 1 )/?(s 2 ) 
- a(s 2 )/?( Sl )] 

The first 3 states form a triplet the total spin is 1. The last is again a singlet 
(like the ground state). 

D. The Slater wave function is the basic anti- 

symmetric function describing N electrons in 
N orbitals 

The previous example is difficult to generalize. In order to develop a way to 
easily represent antisymmetric functions of all types, we consider the 
following 2-electron function, composed of 2 1-electron spin-orbitals: 
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1 

X ¥(X 1 ,X 2 ) = ^=[0i(x 1 )0 2 (x 2 ) - 0i(x 2 )0 2 (x 1 )] 



1 



01 Ol) 01 (* 2 ) 
02 (*l) 02 O2) 



(3.4.1) 



(3.4.2) 



If we choose the orbitals to be orthonormal, (0j|0 y ) = ffy then: 
<V|V> = If¥(x 1 ,x 2 )*^(x 1 ,x 2 )dx 1 dx 2 

= ^ jj [0i(^i)02fe) - 0i(x 2 )0 2 (x 1 )] 2 dx 1 dx 2 

= \\\ [01(X!) 2 2 (X 2 ) 2 + 1 (X 2 ) 2 2 (X 1 ) 2 

+ 20 1 (x 2 )0 2 (x 1 )0 1 (x 2 )0 2 (x 1 )]dx 1 dx 2 = 1 



£?. Without loss of generality, we may assume the 
orbitals of a Slater wave function are 
orthogonal 

But what happens if the orbitals are not orthogonal? Suppose that the orbitals 
were not orthonormal: 

(<Pi\<Pj)=Sij (3.5.1) 

It is then possible to "orthonormalize" them, i.e define two linear 
combinations which are orthonormal. Define: 



fn=^0 



,A 



m ll mn 



(3.5.2) 



and demand: = 8ij> Then: 

fyteto) = Yu A * mi M*j) A *J = ^ SA ^ (3.5.3) 

m,n 

Thus: A^SA = I or AA^ = 5" 1 . Note also that deti4 = (detS) _1/2 . There are 
many solutions to this equation. (For example, if A is a solution then so is 
AqU where U is any unitary matrix.) Each solution will give us a different set 
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of orthonormal orbitals. The Slater wave function made out from these new 
orbitals is: 

detffcfo)} = det{G4 r ) n£ £ (x ; -)} = detA x detfafa)) 

1 (3.5.4) 

= VHiff det{ * i(x ' )} 

Thus, the new wavefunction is the same as the old one, up to multiplication 
by constant! Yet, it is always more convenient to work with normalized 
orbitals, so we can assume the orbitals are orthonormal without any loss of 
generality. This development also shows that given any set of N orbitals from 
which the Slater wave function has been constructed, we can take N linear 
combinations of the orbitals to obtain new orbitals that give the same Slater 
wave function up to a constant factor. 

F.Any antisymmetric function can be expanded as a 
sum of basic Slater (determinantal) functions 

For orthonormal orbitals, the normalization is easy to compute. We write 
explicitly the determinant as: 

N 

det[0! ... 4> N ] (*! x N ) = ^ a-f^N e Y[ cp ik {x k ) (3.6.1) 

t±..J.N k = l 

Where (ii..t w ) is a permutation of the numbers 1 ... iV (there are N\ such 
permutations). Each permutation can be obtained from a series of pair 
swapping operations. For example: (132) is obtained from (123) by switching 
the pair of numbers in position 2 and 3. We write this as: (132) = S 23 (123). 
(2413) is obtained from (1234) by three operations: 

S 23 S 34 S 12 (1234) = S 23 S 34 (2134) = S 23 (2143) = (24 1 3) (3.6.2) 
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If the number of switches is odd the permutation is odd and Pi 1 _j N = 1; if the 
number of switches is even, the permutation is even and P^...^ = 0. The 
normalization of a determinantal wave function composed of orthonormal 
orbitals is: 



| |det[0! .. 



N N 

= £ ( -f^ Ne( _fh.J Ne jY\cP ik (x k )Y\^ jq M dX l- dX N 

N 



(3.6.3) 

ii-ijvA-Jjv fc=i <?=l 



ix-iNh-jN fc=1 
Because of orthonormality the orbital integral (0tj0; k ) is zero unless 4 = j k . 
These integrals appear in products so the product is non-zero only if i k = j k 
for all k = 1 ...N. The only conclusion is, that the two permutations must be 
identical and: 



||det[^ ...Z N ](x 1 ,...,x N )\ 2 dx 1 ...dx N = ^ 1 = JV! (3 64) 

i\...iNH~iN p 



We conclude that the normalization factor of a determinantal wave function of 
orthonormal orbitals is -p= and write: 

1 :det[^ i ...^ N ] = |i 1 ...i w ) (3.6.5) 



Given a set of M > N orthonormal single-electron spin-orbitals (p n (x) = 
(p n (r , s) n = 1,2, ... , M, we can consider the space of all linear combinations of 

all N -particle determinants that can be made. There are ^ ^ = Nl ^'_ N y ways 

to select determinants so this is the dimension of the space. The dimension 

grows factorially with M. A typical antisymmetric wave function can e 

approached by linear combinations of these determinants: 
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W(x 1 ,...,x N ) = -= ^ C k iN det^OO ...<p iN (x N )] ( 3 66 ) 

" {tl tw) 

The sum is over all selections of N integers, where each selection is ordered so 

that i x < i 2 < If the orbitals are orthogonal, the constants C tl iN are 

obtained from: 

C k i N = -^=J det^C^) ...<p iN (x N )]V(x 1 ,...,x N )dx 1 ...dx N (3.6.7) 

G. Determinant expectation values 

In this section we discuss the calculation of expectation values of many- 
electron operators for N electrons within a given Slater wave function 
^sC*!, —,x N ) = ^Ldet^Oq) ...(p N (x N ) ]. We assume the orbitals 0j(x) are 
orthonormal: (0 £ |0y) = Sy. 

i. One-body operators 

Consider an operator d^x) which operates on an electron with spin- 
coordinates x. For N electrons we define the sum of d t for each electron 

N 

= £ d^x n ), (3.6.8) 

n= 1 

Examples: when electrons are in a potential well v(r), the total potential 
energy operator is V = T.n=i v (. r n)> the total kinetic energy is: 



When the system of N electrons is in a given Slater wave function 
^sOt-L, ...,x N ) = ^zdet^!^) ...<p N (x N ) ], then using the notation of (3.6.3), 
we have: 
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(Vs\6\v s ) 

= Jnl det ^^ x ^ -^w(^w)] det^C^) ...<p N (x N )] dx 1 ...dx N 

(3.6.9) 

4zz (-^n^i^^^ 1 ^) 

n=l ii-.ijv fe=l 

Once again, a massive cancellation of terms happens in the first integral. 
Inspection shows that both permutations, i and must be equal otherwise 
there is always an orbital integral for which the integrals (ft fe hy k ) is zero. 
When the permutations are identical we have: 

N N 
n=l i x ...ijv n=\ 

As an example, let us take the electron density n(r) = £n=i S(r — r n ). Thus: 

N 

n{r) = Y\<Pn(X,s)\ 2 (3.6.11) 

n=l 

Conclusion: The matrix element of a 1-particle operator is the sum of its 
single-electron matrix elements. 

ii. Two-body operators 

Consider an operator d 12 (x 1 ,x 2 ) which operates on two electrons with spin- 
coordinates x 1 and x 2 . For N electrons we define the sum of o 12 on all pairs of 
electrons 

N N 

= 2^ (*n-*m) = 2 ^ o 12 {x n ,x m ), (3.6.12) 
n<m= 1 n*m= 1 

Examples: The 2-electron interaction potential is Ui 2 (Ti, ^2) = r~T — r The total 



rv 1 ■ 

2 ' 



interaction energy operator is: U = -Ef=iEf=i lt i2( r t< r ;)- 
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We compute the expectation value 

= ^ j det [0i(*i) -0n(*n)] U det^C^) ...(p N (x N )] dx 1 ...dx h 

N N 



N 

1 1 



(3.6.13) 



JV!2 . 

n,m=l i 1 ...i N k=tn 
n=tm j....j N k±m 



Where we used the notation: 

(0i0; K 2 \<Pi ,( Pj') = |0i(x)0 ; (x')ui 2 (r,r')0i'W0y'(^')d^' (3.6.14) 
The following symmetry properties hold from the above definition: 



(0i0y|tll2l0i'0/) = (0;'0;K 2 |0 £ </y) = (0i0/K 2 |0 £ '0;) 
= (0£'0 ; '|"i 2 l0i 0;) = etc 



(3.6.15) 



For a pair of permutations to contribute to the integral in Eq. (3.6.13), the 
permutations must either be identical or involve the permutation of a single 
pair of orbitals. Thus: 

1 N 

(V S |^|V S ) = - ^ [( <p n <Pm\u 12 \<Pn<Pm) ~ ( <Pn<Pm\u 12 \<Pm<Pn)] (3.6.16) 



n,m=l 



IV. The Hartree-Fock Theory 

A. The ffartree-Foak Energy and Equations 

The variational principle says that the lowest expectation value of the 
electronic Hamiltonian attained by the ground-state. This Hamiltonian, in the 
non-relativistic approximation, for N electrons is given by: 
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N ,2 N 



n=l n*m=l 



(4.1.1) 



Where the "one body" operator is 



^ = Z^ = Z _ 2^: V ^ + u(r?n) (4 - L2) 

m=l m=l 

And the two-body Coulomb repulsion operator is: 

N 

I/ = T > i r 4.1.3) 

2 La \r n -r m \ 

n*m=l 

Given a family of wave functions we can choose the "best" of them by finding 
that which minimizes the expectation value of H. For the Slater wave 
functions W s = |0 X ... (p N \, the energy to be minimized is: 

E[V S ] = (V s \8\Vs) = EhfWi -0n] (4-1.4) 
Because we want the orbitals to be orthonormal, we write a Lagrangian: 

N 

L[V S ] = ...0 W ] = £ HF [0! ...0 N ] - £ IHj[(4>i\4>]) - (4-1.5) 

The minimum is attained by: 



5L Hf [0! ...0 W ] 5E HF [0! ...0 W ] , \. ,^ 



(4.1.6) 



;'=i 

Let us try a solution with 



i.e.: 



= CifaQc) (4.1.8) 



8</>i(x) 
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If we can solve this equation, and if the solutions are naturally orthogonal, 
then we have obtained the necessary conditions for a minimum. Now, we 
only need to estimate the left hand side of this equation. From the previous 
work, we know: 



(4.1.9) 



where 



N 



(V S \h\V S ) = h[01 -0 W J = ^(0mN0m) 



(4.1.10) 



771=1 



Is the "one body energy" defined the as the sum of kinetic energy and 
"external" energy (i.e. energy due to the frozen nuclei): 

h 2 



h t = - - — V 2 + v(r) 

Z[J. e 

Furthermore, the 2-electron energy can be written as: 



(4.1.11) 



(V s \0\V s ) = e - ]T (</> n r 



71,771 = 1 
77*771 



0770777 - 0770, 



12 



0777077 



(4.1.12) 



=y[0i-0 w ] + ^[0i-0 w ] 

Where the direct or Hartree energy is: 



N 



7[01 -0W] = y V (077071 

and the energy exchange is 



77,777 = 1 









(0770777 




07707n) 



(4.1.13) 



77,777 = 1 





1 




(0770777 


ri2 


0777077^ 



(4.1.14) 
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The direct energy is numerically equal to the Hartree energy E H [n] which is a 
functional of the Slater wave function electron density n(r) = (^InCr)! 1 !^) = 

Zm=il0mfcs)| 2 : 

BM^f^™, (4.1.15, 

The exchange energy is numerically equal to the "exchange energy" E x [p] 
which is a functional of the density matrix defined by 

N 

P {X,X') = 0m(*)0mOO (4-1-16) 



m=l 



As: 



\p(x,x')\ 2 



dxdx' (4.1.17) 



\r — r'\ 

Notice that: j p{x, x)ds = n(r). The density matrix is idem-potent: 

jp(x,x"Mx",x')dV = p(x,x') (4.1.18) 

This result shows that the DM is a projection operator, projecting onto the 
space of orbitals which defines the Slater wave function. 

With the direct and exchange energies we also define their functional 
derivatives: 
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5] 



1 Y 8 



50*00 2^50*00 





e 2 




{$n<Pm 


ri2 





= yl^JJ |r t -r 2 | 

e 2 v / r 0m(^2) 2 



+ 28 im (f) m {x) J _ ^ I 

= 2e2 (/ |r-r 2 | d3rz )^ f ^ S 2v H (r)<p i (x) 



= 2v H <pi(x) 

Where in the last line we defined the Hartree potential: 



(4.1.19) 



v H (r) = e- 



ZjJ |r-r'| 



J |r-r'| 



m=l 



(4.1.20) 



Then the direct energy can also be written as a functional of the orbitals: 

N 

;[0i-0 W ] = 2 X <0ml ^ l0m> = 2 J ^(r)n(r)d 3 r = -(^ s |K H |^ s ) (4.1.21) 

m=l 

where = £n=i u H( r n) is the total Hartree potential. A similar treatment 
exists for the exchange energy functional derivative: 

SK e 2 Y s 





1 




(0n0m 


^12 


0m0n) 



-y^[25 in m (^) J 



+ 28 im (p n (x) j 
= -2e 



\r-r 2 \ 

<Pn(x 2 )<p m (x 2 ) 



(4.1.22) 



\r-r 2 \ 



■dx 7 — 



2 y f 0nfe)0i(^ 2 ) 

2J |r-r 2 | 



dr 2 (p n (x) = 2K(pi(x) 
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Where the last equality is based on a definition of a one-particle exchange 
operator: 

N 

KupCx) = -e z 2^ <t>m (*) I . _ r ,, — dx' 

™=1 " (4.1.23) 

c p{x,x')\l){x') , 

= -e z — = — dx 

J \r-r'\ 

Then the exchange energy is written as a functional of the orbitals: 

N 

K[<t>x = \Y J ^n%\^>n) = \ ^s\Qh\Vs) (4-1.24) 

n=l 

where Q H = Y,n=i Kn is the total exchange operator. The other functional 
derivatives needed are: 



Sep 

and 

o 



8 C 

— (P n (x')V 2 (P n (x')d 3 x' = 2V 2 <p n (x)5 in (4.1.25) 
iW J 



U f 

^~r^ <t>n(x') 2 v{r')dx' = 2v(j)<p n {x)8 in (4.1.26) 
Thus: 

— ^— f cf> n (x')h(r')(i> n (x')dx' = 2hcp n (x)S in (4.1.27) 

Plugging all these terms into Eq. (4.1.7), we obtain the Hartree-Fock 
equations: 

F<f> t (x) = erfiix) (4.1.28) 

where: 

F = h + v H (r) + K x (4.1.29) 
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Equations (4.1.28) seem very much like 1-electron eigenvalue equations of the 
Schrodinger equation, except that instead of a regular Hamiltonian, we have a 
Fock operator including the non-local exchange. We showed in the exercise 
that K is Hermitean and thus so is F, i.e. (0|f|i/;) = (i/;|f|0) . We can thus 
choose the orbital solutions of the eigenvalue equation (4.1.28) orthonormal. 
This shows that the choice Eq. (4.1.7) is indeed acceptable. 

Because v H (r) and K themselves depend on m , the Hartree-Fock equations 
are fundamentally different from the Schrodinger Equation: they are nonlinear 
equations. 

Now that the sum of orbital energies is: 

N N N 

2^e m = ^<0 m |F|0 m ) = ^[((pmlh^cpm) + ((p m \v H +K 1 \4> m )] (4130) 

771 = 1 777=1 777 = 1 V " " / 

= h[<p lt ... , W ] + 2/[0 1 , ... , W ] + 2tf [0i, - , W ] 

This shows that the orbital sum is not equal to the energy of the wave 
function, since it involves double counting of the direct and exchange 
energies. The HF energy is thus: 

N 

E HF [(pi <t> N ] = Yu £m ~ Wi 0w] + K[<pi 0w]) (4 - L3i) 

777=1 

B . Restricted closed-shell Hartree-Fock 

For molecules with even number N of electrons in a spin-singlet state, we can 

impose the following structure on the Slater wave function. We can assume 

that the 2N spin-orbitals come in pairs: 2 ;-i(X) = 4>j(v)a(o)) and (p 2 j(x) — 

xpj(r)P(a)). Thus each orbital ipj, j = 1, ...,N/2 is "doubly occupied" by 

electrons of both spins. By imposing this constraint we obtain the "restricted" 

Hartree-Fock ground state. It will sometime be of higher energy than the fully 

Electron Density Functional Theory Page 95 

© Roi Baer 



unrestricted case. However, the wave function has a well-defined spin which 
may be advantageous in some applications. 

The Restricted Hartree-Fock (RHF) energy in the closed shell case remain 
essentially the same except for counting business. We can formulate all 
expressions using only the spatial orbitals. Indeed, the RHF energy is given 
by: 

Erhf[*Pi> ■■■'V'w/2] 

^ 2 _ (4.2.1) 
= 2 ^{\p m \h\\p m ) + AJ^, ... ,\p N /2 ] +2K[x(j 1 ,...,x(j n/2 ] 

m=l 

Where ]\ip lt —.ipN/z] an d K.\\p lt ip N / 2 ] are the orbital functionals defined in 
Eqs. (4.1.13) and (4.1.14) respectively. The reason for multiplying the one 
body part by two is evident: each orbitals is double occupied so has double 
contribution. The direct part is multiplied by 4 since the density is multiplied 
by two and the direct part depends on the density multiplied by itself. Finally, 
the exchange part is multiplied by 2 and not 4 since only a — a and /? — /? 
contribute, while a — /? and (3 — a do not (so only half the contribution of 
direct). 

The RHF equations are then: 

( h 2 . \ 
72 



V 2 + v(r) + v H (r) xp m (r) + Krf> m (x) = e m xp m (r) (4.2.2) 



V 2jU 
Where: 

N/2 



v H (r) = 2e 2 \ f 'f ™ (r)|2 d 3 r' (RHF) (4.2.3) 
Z_i J \r — r'\ 

m=l 



and 
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N/2 

K^{r) = -e- £ tf»Cr) / dr' (RHF) (4.2.4) 

m=l 

Note that the Hartree interaction is between each electron and all other 
electrons regardless of their spin while the exchange interaction involves each 
electrons with all other electrons of the same spin. 

Example: The H2 molecule 

We apply the RHF theory for the H2, having a pair of electrons. The 2-electron 
wave function includes just one spatial orbital populated by spin-paired 
electrons: 



0(x 1; x 2 ) 



^(rj xp(r 2 ) 



iKr!>Kr 2 )[a(l)/?(2) - a(2)/?(l)] (4.2.5) 



Since there is only one orbital the exchange K and direct J involve just one and 
the same integral. Thus the RHF energy is in this case: 
E RHF [ip] = 2(ip m \h\if} m ) + 2][\p] and the RHF equation is 

V 2 + v(r) +-v H (r) \tp(r) = ei/j(r) (4.2.6) 



2fi 



where v(r) = —— - — : — — - — r and v H (r) = 2e 2 f ^-^-d 3 r'. The effect of 

v J \r-R A \ \r-R B \ HK J 3 \r-r'\ 

exchange here is to annihilate the Coulomb repulsion of the a (/?) electron 
with itself, leaving only the interaction of the a — (3 electrons. 

This RHF approach works nicely for the case that the distance between the 
nuclei \R A — R B \ is close to the typical bond length of H 2 (which is close to 
1.4a ) . The energy of the molecule at this configuration can be calculated 
numerically and results in E RHF = —1.134E h . Compared to the energy of 2 H 
atoms (— lE h ) this results indicates that the atomization energy of H 2 is 
0.134E/J = 3.65eV. The atomization energy based on expremintal results is 
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about 4.75eV. This shows that the RHF approximation does not give high 
quality atomization energies, since a deficit of 1 eV is very substantial in 
Chemical energy terms. 

A more severe problem arises when we place the two nuclei far from each 
other. We expect the resulting energy and wave function to resemble that of 
two separated H atoms. I.e the exact wave-function should approach 

something like W = ^(Ti ~ Ra)$is(T2 -Rb)\~ |<Ks0 2 - Ra)^u(Ti ~ 
i? B )|.This form however is not supported by the RHF ansatz of Eq. (4.2.5). 
Indeed, if we think of the solution of the RHF equation as being 
approximately given by xp ~ xp A (r) + ip B (r), where ip x (r) = xp(r — R x ), 
X — A,B, then the RHF wave function is: 



The first term is an ionic term, where both electrons are on the same atom 
(either A or B) while the second term places one electron on each atom - 
"neutral" term. The problem of RHF theory is the ionic term. It may be 
important when the atoms are close but it should go to zero when they are far. 

C. Atomic Orbitals and Gaussian Basis sets 

Where do we get "good" basis functions? What is "good"? 

We want a small basis that can still describe the electrons. On natural source 
are the atomic orbitals of the atoms. These are of the form resembling 
exponentials time polynomials. Thus, one choice is: 




= [|Wri)«r 2 )| + hh,(fi)^(r 2 )|] 
+ [|Vu(*i)$B(r 2 )| + |Vs(*i)^(r 2 )|] 



(4.2.7) 



xL<r) = Q&r - R A ) 



(4.2.8) 
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Where: 



% m (r)=r l e-< r Y lrn (e,4>) 



(4.2.9) 



Where Y im are the spherical harmonics. One can also take appropriate 
combinations of these functions to make them all real. These functions have a 
desired analytical property: their Z's derivative exhibits a cusp of the correct 
order and structure at r = R A . There exist analytical formula for doing the 
overlap and "one-body" integrals. But there are no convenient formulae for 
the 2-body integrals, although some progress was made in recent years (see 
articles by Handy). 

A more convenient, although less natural choice (no cusp). Is the use of 
Gaussian functions, for example: 



Where a ifi are called "contraction coefficients. These are chosen so that 



rapid algorithms were published allowing extremely fast 2-electron integrals. 

D. Varlational-Algebraic approach Hartree-Fock 

We have seen that the Hartree-Fock equations can be derived by searching for 
that the most general Slater wave function that minimizes the Hartree-Fock 
functional. However implementing a solution to such equations is usually 
very difficult, if not impossible in practice. A more practical approach, that 
keeps the spirit of the Hartree_Fock approach was developed by Roothan and 
Hall. In this approach we find the optimal Slater wave function of orbitals 
constrained to lie in a finite dimensional vector space spanned by basis 




(4.2.10) 




) resembles e ^ r . With Gaussian functions very effective and 
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, M 


M i 


Snm = 1 


^ ' Xa^an 


^ Xa' C a'mj 









functions , usually called atomic orbitals (although, up to a point, we need not 
assume this) Xa(. x )> ° — 1, Thus, a set of N molecular orbitals (p n (x) 

(n = 1, ... N) in a determinant ^[Cjof this form must all be of the following 
form: 

M 

4>n(x) =Y 4 XMC an (4.3.1) 

(7=1 

The C an coeficients form anMxW matrix, called the MO coefficient matrix 
for the determinant. Note that for this to make sense we must demand M > N. 
The HF energy functional now becomes a function of these coefficients C an . 
The constrained that the MO's are orthonormal, {(p n \<p m ) = S nm , becomes: 

(C T SC) nm (4.3.2) 

Where we use matrix algebra notation and the M x M matrix S is defined by: 

Soo> =iXo\Xa>) (4.3.3) 

Thus, the orthonormality condition is: 

C T SC = I N (4.3.4) 

Where I N is the N x N unit matrix. Let us now derive an expression for the 
expectation value of a one-body operator in a Slater wave function of these 
MO's, by Eq. (3.6.10): 

N 

(^ 5 |0|¥ s ) = ^(0 n |O 1 |0 n ) = Tr[C T 0C] (4.3.5) 

n=l 

Where is the MxM matrix in the AO basis: 

Ooo' ={Xa\6 1 \Xo<) (4-3.6) 

It is customary to define the MxM density matrix: 
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P = CC T (4.3.7) 

And with it we can write: 

(V 5 |0|V s ) = Tr[C T OC] = Tr[PO] (4.3.8) 

In the last step we used the fact that the trace of the product of two matrices is 
invariant to their order of multiplication: 

M M N N M 

Tr[AB] = ^045)^ = ^ ^ A an B na = ^ ^ B na A an = Tr[BA] (4.3.9) 

cr=l cr=l n=l n=l cr=l 

Notice that the DM has the generalized idempotency property: 

PSP = CC T SCC T = CC T =P (4.3.10) 

One can see that P is a symmetric matrix. Furthermore, one can see that it is 
positive semi-definite, i.e. for any vector: v T Pv = v T CC T v = (C T v) T (C T v) > 0. 
Furthermore: 

Tr[PS] = Tr[CC T S] = Tr[C T SC] = N 

This last step is a result of Eq. (4.3.4). Finally the 2-body operator, by Eq. 
(4.1.12), we need the direct and exchange. We use: 

N 

[kl\mn] = ^ C T ka C all [oo'\TT'}Cl T C T , n , (4.3.11) 



Then: 



e 2 



(4.3.12) 



e 2 



K[C] = -Y C n*C a >m[°°'\ TT '] C LC T , n , 
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Where we use the convention that all repeated indices are summed over: latin 
indices are summed between 1 and N and greek indices between 1 and M. The 
summations on n, m can be done first and we obtain: 



e 2 



e 2 



(4.3.13) 



K[C] = - — P a A™'\™'\Po'i 



The last expression can be reindexed (assuming P is a symmetric matrix) as: 



e 2 



K[C] = --P a , (J [GT'\TG']P W (4.3.14) 



Thus: 



e 2 



J[C]+K[C] = -P a > a {[aa'\rr'] - [ot'\to'\)P tt , (4.3.15) 
Finally, defining: 

Vij = V aa , iTT , = [aa 1 \tt'] - [ar'\Ta'] = [aa 1 \tt'] - [ar'la'r] (4.3.16) 

Using the double indexing / = (erer") and / = (o'o'"). Note that V t] = V ]l as 
can be seen from: 

V u = V aa , iTT , = [aa' \tt'] - [ar'la'r] = [rr'\aa'] - [ra'\r'a] = V TT , aa , 



= V ]I 



(4.3.17) 
We may thus write: 



e 2 



J[C] + K[C] = —P l VP (4.3.18) 



Where now, we consider P not as a M x M matrix, but as a column vector of 
M 2 elements. P t is the corresponding row- vector. Similarly, V is not a M x M X 
M x M tensor but as M 2 x M 2 matrix. 
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The HF energy is thus compactly written as: 

1 

E HF [P] = P f /i + -ptVP (4.3.19) 

Note, that our unknown variable is now the DM P. We want to minimize E HF 
with respect to P, however, we need to impose 2 types of constraints. First, we 
need to specify thesubject to the constraints: 

Tr[PS] = N 

(4.3.20) 

G = PSPS -PS = 

In order to minimize the energy we introduce the Lagrangian: 

L[P] = E HF [P] - n(P'S - JV) - A t G (4.3.21) 

The number n and the M x M matrix A are Lagrange multipliers. The 
algebraic Hartree-Fock equations are now = 0. In order to obtain working 

"Par ' 

expressions we derive: 

SE 1 

~spy = h ' + 2 ( p M' + v v p *) = h ' + Vl > p > (43 - 22) 

Where again, we use the convention that when a super index appears twice 
we sum over it. This can be written more compactly as: 



SE HF 



h + VP = F (4.3.23) 



SP 

This gradient is what we call the "Fock matrix" F. In our present notation F is 
a M 2 vector Fj with index /. But soon we will consider it as a matrix with two 
indices F aa >. The constraints can easily be derived in a similar way, leading to 
the following Lagrangian gradient: 

— = F-[iS- (SPSA + SAPS - SA) (4.3.24) 

It is easy to convinve one's self that F is always a symmetric matrix, for any P. 
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The condition for minimum is: 

F-fiS- (SPSA + SAPS - SA) = (4.3.25) 
Multplying by SP from the left we find: 

SPF - [iSPS - SPSAPS = (4.3.26) 
Multplying by PS from the right we find: 

FPS - fiSPS - SPSAPS = (4.3.27) 
Subtracting, we obtain: 

SPF -FPS = (4.3.28) 
The set of equations that need to be solved simultaneously is: 
SPF - FPS = F = h + VP 

(4.3.29) 

PSP = P Tr[PS] = N 

One practical way of doing this is to go back to the matrix C. In terms of these, 
the equations become: 

SCC T F - FCC T S = F = h + VCC T C T SC = I N (4.3.30) 

These equations can all be met if we demand that: 

FC = SCE F = h + VCC T (4.3.31) 

Where E is a M x M diagonal matrix. Indeed, from this equation we also have, 
from the symmetry of F and 5: C T F — EC T S. Left-multiplying by SC we find 
SCC T F = SCEC T S and using the first equation in (4.3.31) on the right hand 
side we obtain the first equation in (4.3.30). Furthermore, multiplying the first 
equation in (4.3.31) by C T we find: C T FC = C T SCE. On the left we replace C T F 
by EC T S and obtain: 

[E, C T SC] = 



Electron Density Functional Theory 
© Roi Baer 



Page 104 



We find that C T SC is commutative with a diagoinal matrix. If no two elements 
on the diagonal of E are equal then C T SC is diagonal. We know that the 
diagonal entries must be positive since S is positive definite. Furthermore, we 
can choose the norm of the columns of C so that all diagonal elements of C T SC 
are equal to 1. In this case then C T SC = I. When there are several elements on 
the diagonal of E which are exactly equal, then one can take linear 
combinations of the corresponding columns of the C-matrix, without 
disturbing their eigenstatishness. Once can then always create a situation 
which again allows C T SC = I . We thus find that the procedure of finding the 
generalized eigenstates and eigenvalues of F is indeed a procedure for finding 
the minimum. 

Thus Eq. (4.3.31) is the algebraic Hartree Fock equation. In actual calculations, 
it is very common that programs solve self consistently the algebraic HF 
equation. This procedure is appropriate for small to medium sized systems. 
But for larger system it may be beneficial to directly the minimize of the 
Lagrangian, using the gradient in Eq. (4.3.24). Of course, iterations are still 
needed because a search must be made for the Lagrange multiplier A. 

E. The Algebraic Density Matrix and Charge 
Analysis 

We have seen that the density matrix is defined by the relation P = CC T where 
C an is the coefficient of the AO Xa( x ) i n the expansion of the MO <p n (x). The 
relation of P the the real space density matrix defined in XXX is (we use the 
convention that repeated roman indices are summed from 1 to W and 
repeated greek indices are summed from 1 to M): 

P(X,X') = n (*)0 n (*') = C an X a {x)Xp(x')C pn = PapXaWXptx') 
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Thus we see that P a g determines p. Hence the common name for these two 
quantities. Note also that P aB determines the density n(r) since that is easily 
obtained by place x = x' in the DM and integrating over spin: 

n(r) = hm + P aB S aB — — ■— - = q a n a (r) + ^ q a pn aB (r) 

Where: 

x a (r) 2 



n a (r) = 



S 

°aa 



Are the "atomic electron number densities" (each integrates to 1) and 

. . Xa(r)xp(r) 

n aP (r) = 

Is the bond electron number density (again, integrating to 1, or zero). The 
atomic charges are then q a = P aa S aa and the bond charges areq ap = P a pS ap . 

This form of charge analysis is very popular and allows to obtain" intuitive" 
pictures for the charge distribution in the molecule. While useful to many the 
user should be warned that this analysis is "basis-form" dependent. What we 
mean by this is that if we take different linear combitaions of the same set of 
basis functions (i.e. we stay in the same Hilbert space), our charge analysis 
will yield totally different results. This is because when we take linear 
combinations: x' = then the density matrix changes by P' = T T PT. Thus in 
general the charges on each atom can change by this procedure. 

Note that when one integrates over r , on the left hand one gets N. On the 
right hand the first gives (assuming the basis functions are normalized P aa 
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F . Solving the Hartree-Fock Equations 

A plausible "algorithm" for solving the Hartree-Fock equations is as follows: 

1. Guess <pm(r), m = 1, ... , N e . 

2. Build v H (r) (Eq.(4.1.20)), K (Eq.(4.1.23)) thus determining F (eq. 
(4.1.29)). 

3. Solve the eigenvalue equations (Eq. (4.1.28) to get a new set of orbitals 
corresponding to the lowest energy orbitals. 

4. Redo from step 2 using the new orbitals, until you converge - i.e. until 
the orbitals change no more. 

While this algorithm seems reasonable, in practice it rarely converges. There 

are several ways to make an algorithm "practical", 
i. Direct inversion in iterative space (DIIS) 

This method, devised by Pulay (P. Pulay, Chem. Phys. Lett. 73, 393 (1980)) is 
designed to speed up the convergence. Suppose the iterative process has 
produced M iterants v m m = 1, ... , M (Fockians, density matrices or sets of N e 
orbitals). We can define residuals by: 

8v m = v m - v m _ x (4.4.1) 

We want to produce a new iterant by interpolation: 

M 

v = ^ w m v m (4.4.2) 

m=l 

where w m are the weights and they sum to unity: 

M 

£ w m = 1 (4-4.3) 

771=1 

These weights are obtained by minimizing the residual, assuming linearity: 

M 

= ^ w m Sv m (4.4.4) 



M 

Sv 

771 = 1 
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The function to be minimized is: 

M 



J\w] = —Sv T Sv — X 
2 



I 

m=l 
M 



w m -l 



= \1 



w rrfim,m' W m' ~ ^ 



Where: 



m,m'=l 



B m,m' = 5 ^m S ^m' 



M 

I 

m=l 



w m -l 



(4.4.5) 



(4.4.6) 



Differentiating with respect to w k gives: 

dj 







dw k 



Bkm w m ~ ^ w k 



(4.4.7) 



The solution of these equations, together with the constraints Eq. (4.4.3) gives: 





B n ■ 


B IM 


-r 




fw 1 " 






B 21 


B 22 ■ 


■ ■ B 2M 


-i 




w 2 









B M2 


B MM 


-i 




w M 







1 


1 ■ 


■■ 1 







A 




1 

















(4.4.8) 



The solution of this equation gives the desired weights. The use of this 
algorithm can be done in the following way: 

1. Get new v as output from the iterative procedure. Add it to the list i.e. 
designate it as v M 

2. Find weights from which get interpolant v' = Xm=i w m v m . 

3. Use v' as input to the iterative procedure and redo from 1. 

ii. Direct Minimization 

Sometimes the DIIS procedure is not effective and other methods are tried. 
One of the most useful methods is to use numerical minimization techniques, 
such as the conjugate gradients algorithm to directly minimize the energy of 
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the Slater wave function under the constraints. This methods is especially 
useful when the molecule being studied is very large. Special tricks are used 
to formulate the minimization problem is an "unconstraint minimization" 
(see for example, Nunes et al, Phys. Rev. B 17611, 50 (1994)). 

G. Performance of the Hartree-Fock approximation 

We examine the performance of Hartree-Fock approximation on, for example 
Formaldehyde. There are 2 sources of error. One is in the application, since we 
use finite basis sets. Then there is the intrinsic error. 

In the table below, we see the prediction of various properties of 
formaldehyde, calculated with increasing quality of basis set and compared to 
experiment. 



Basis 




R(CO)A 


R(CH)A 


A(O-C-H) 


Energy(au) 


sto-3g 




1.2169 


1.1014 


122.73 


-112.35435 


3-21g 




1.2071 


1.0833 


122.51 


-113.22182 


sto-6g 




1.2163 


1.0981 


122.61 


-113.44078 


6-3 lg 




1.2103 


1.0816 


121.69 


-113.80837 


D95 




1.2170 


1.0843 


121.57 


-113.83071 


D95v* 




1.1887 


1.0935 


121.96 


-113.89173 


6-31 lg** 




1.1787 


1.0949 


122.09 


-113.89915 


6-31 l++g** 


1.1797 


1.0943 


121.97 


-113.90287 


apvtz 




1.1786 


1.0927 


121.94 


-113.91534 


experimental 


1.210 


1.1020 


121.1 





We see Hartree-Fock converges when basis set quality increases. However the 
converged quantity deviates somewhat from experimental values. 

This deviance exists because Hartree-Fock theory is only an approximation. 
What it assumes is that the electrons act as if they are independent particles 
(since it imposes a single determinant). The real ground-state is composed 
from a huge series of determinants. The "independent" particles interact with 
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the mean field of all other particles, while in essence each electrons has to 
interact with each other electron, trying as much as possible to avoid it, 
without paying too much in kinetic or electron-nuclear potential energy 

H. Beyond Hartree-Fock 

The Hartree-Fock method is very successful, since it typically accounts for 
over 99% of the electronic energy of molecules. Yet it is not accurate enough 
for most applications in chemistry. The reason is that most quantities of 
chemical significance are energy differences - not absolute energies. When 
differences are considered the errors in the Hartree-Fock approach are not 
small. 

One way to improve the situation is to approximate the groundstate wave 
function by a series of determinants: 



where (p n are an infinite orthonormal set of orbitals. Such an expansion can 
always be made, with any such set. We can thus take the orbitals produced by 
the Hartree-Fock process. This has the added nicety that the first determinant 
is already a good approximation to the ground-state. 

In this case we can classify the determinants in the following way. We divide 
the orbitals into two sets: one is the set of N e HF orbitals, called the occupied 
orbitals, the rest of orbitals named virtuals. We then classify the determinants 
by the number of occupieds missing. Thus we speak of all single substitutions, 
double substitutions etc. 
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(4.6.1) 



A commonly employed method is the configuration-interaction (CI) method. 
In a CI one takes a number of determinants D a a = 1, ...,M and uses them to 
minimize the energy: 

M 

E[c] = £ c m {D m \H\D n )c n (4.6.2) 

m,n=l 

under the constraint of normalization: £m,n=i c m (D m \D n )c n = 1. One common 
way of choosing the determinants that go into this expansion is by collecting 
all single, double, triple etc excitations. A determinant is singly excited if 
when compared to the HF determinant it has one occupied molecular orbital 
replaced by some virtual orbital. Virtual orbitals are excited eigenfunctions of 
the Fock operator. One can show that the singles alone do not allow a 
reduction of energy. However, singles and doubles give sometimes good 
results. Such a method is called singles-doubles CI (SDCI). One problem with 
this theory is that it is not "size consistent". For example, calculating the 
energy of 2 distant Helium atoms will not give the twice the energy of one 
Helium atom under the same order of theory. 

Another approach is many-body perturbation theory, called Moller-Plesser 
theory. In this approach, one writes the many-body Hamiltonian as: 



M 




(4.6.3) 



n=l 



Where W = H — En=i Fn is considered a "small" perturbation. This quantity is 
not small enough and high order MP theory does not converge. However, 
second order MP theory, called MP2, is sometimes a useful approach. It is size 
consistent. However, it relies heavily on the quality of the Hartree-Fock 
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solution: Hartree-Fock must be a good reference n which to base a 
perturbation theory. 

All wave-function methods beyond Hartree-Fock theory become quickly very 
expensive as system size grows. In fact, the numerical cost of good methods 
typically scales as 0(iV 7 ) for the coupled cluster method, which is a size- 
consistent non-variational method (variational methods are derived from the 
variational theorem), not discussed above. Thus, every enlargement of 
number of electrons by a factor of two makes the calculation a factor 100 more 
expensive! 

V. Advanced topics in Hartree-Fock 
theory 

In this chapter we will continue our study of the Hartree-Fock approximation, 
and look into some of the formal issues, like stability, excitations and 
ionization and generalizations like fractional occupations numbers. We will 
then discuss the homogeneous electron gas in the Hartree-Fock 
approximation and show that it breaks down when treating this system. 

A. Low-lying excitations and the stability of the 
Hartree-Fock ground state 

i. Cl-Singles and Brillouin's theorem 

For a given system, one can think of the simplest excited states as linear 
combinations involving low lying energy determinants. For example, all 
determinants where one occupied orbital in the HF determinant is replaced 
by an orbitals which is a "virtual" eigenstate of the Fockian (we call those N 
eigenstates of the Fockian which are part of the HF determinant "occupied" 
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and those which are not we term "virtual"). In the algebraic form of HF 
theory, if there are N electrons and the size of the basis is M > N then there are 
N ways to choose the occupied orbital to be replaced and (M — iV) ways to 
choose the virtual orbital and so there are iV(M — iV) such "singly excited 
determinant". Together with the HF determinants we can form a N(M — N) + 
1 dimensional determinant space and diagonalize the exact Hamiltonian in it. 
This will give an approximate description of the low lying excited states of the 
system. In fact, there is no need to include the HF determinant itself in this 
scheme since Brillouin showed that for an singly excited determinant (this 
notation is for a determinant that is obtained from the HF determinant by 
replacing the occupied a orbital by the virtual n orbital): 

(¥£|//|¥ ) = (5.1.1) 

Thus, the HF determinant is decoupled from the singly excited determinants 
and one can just diagonalize the Hamiltonian in the singly excted space. This 
approach is often called "configuration interaction - singles" (CIS) and it is a 
standard method for calculating excitation energies in HF theory. In fact, one 
cannot expect CIS to give a good approximation for the excite states, since the 
"real" excited states are intricate linear combinations of determinants with 
multi-excited electrons. On the other hand, our approximate ground state 
wave function is also just a single (Hartree-Fock) determinant and it too 
misses all this essential mixing with doubly and triply etc excitations. One can 
hope however for a mutual cancellation of errors. Indeed, there are many 
examples where the CIS method gives quite respectable excitation energies, 
even when the wave functions are of questionable quality. 
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ii. Hartree-Fock Stability 

We now consider the question of stability of the Hartree-Fock solution. 
Basically the issue is this: how do we know that we have produced the truly 
lowest energy by occupying the N orbitals with the lowest orbital energy? 
Maybe if we used a different ansatz we could have produced a lower energy. 
In other words, how do make sure that all singly excited determinants are 
higher in energy than the HF determinant. 

Suppose we have set up the Hartree-Fock equations and solved them to 
obtain a Fockian, a set of orbitals and orbital energies, out of which N e are 
"occupied". Suppose the Hartree-Fock determinantal wave function is 
= —> 4>a> — > 4>n I- Let us consider the determinantal wave function 
— I0i' ■■■ ' 07i' ■■■ ' 0w e I obtained from the Hartree-Fock function by replacing 
an occupied spin-orbital (p a by the unoccupied spin-orbital (p n . This can be 
viewed as an electron excitation process: a hole is made in cp a and an electron 
is formed in cp n . The excitation energy for this excited state is the difference 
between the expectation values: 

^ = (¥ a "|H|^)-(%|H|¥ ) 

= {<Pn\h\<Pn) - (0a|^l|0a) 

+ ^Vm(l - 8 ma )({nm\nm) - (nm\mn)) (5.1.2) 

m 

— ^ f m ((am\am) — (am\ma)) 

m 

Here, we introduced, for convenience the orbital occupation f m = 1 for 
m < N e and f m = otherwise. Rewriting: 



Electron Density Functional Theory 
© Roi Baer 



Page 114 



e n = 



((p n \hi\<Pn) + ^fm(.(nm\nm) - (nm\mn)) 

m 

{(pa\hi\(p a ) + ^fm((am\am) - (am\ma)) 



(5.1.3) 



— (Jjta\na) — (na\an)) 

Using the definition of orbital energies in (4.1.30) we can thus write this 
"electron hole" excitation energy as: 



where: 



€ n — € — £ — A 
c a c n c a "an 



A an = (an\an) — (an\na). 



(5.1.4) 



(5.1.5) 



One way to understand A an is as an over-counting term. We thus see that 
excitation energies as calculated as the difference between excited state 
energies in Hartree-Fock theory are not simply the differences of the Hartree- 
Fock orbital energies. They must actually be corrected for over-counting by 
subtracting a quantity A an/ which we show henceforth to be manifestly 
positive. One conclusion is that orbital energy differences in Hartree-Fock 
form an upper-bound to the excitation energies, as determined from the 
Hartree-Fock single excited determinants. We may think of — A an as the 
Coulomb energy of attraction between the excited electron and the hole it leaves 
behind. Indeed A an is composed of the electrostatic interaction of the electron 



(p n and the hole <p a , (0 a n — 0a0n) corrected by a corresponding exchange 



\0a0n — 0n 0a /term. 



To prove A an is positive, note that it is symmetrical: k an — A na . Thus: 
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[A an + A na ] 



1 




e 2 




[(pa<Pn 


Tx2 


2 







e 2 




+ {$n<Pa 




<Pn<Pa 








1/ 






= 2Y a< ^ 


n ~ 


<Pn<Pa 



4>a ( Pn ~ <Pn<Paj 

0a0n - 0n0a i 



(5.1.6) 



'12 



This shows that 



A an =y// 



e 2 ff (0 a (l)0 n (2)-0 n (2)0 n (l)) (3 w3 
\r 2 -r x \ 



d 3 r!d 3 r 2 > 



(5.1.7) 



Stability is obtained when > 0, or, in other words when e n — e a > Since 
A" > 0, we find that a necessary condition for stability is: 

e n - e a > (5.1.8) 

Clearly this might not be sufficient, however, it is necessary . That is: for 
stability to be possible, the orbital energy of all occupied orbitals must be 
lower than that of all the unoccupied orbitals otherwise there will definitely 
definitely be "singly excited determinants" with lower energy. 

B.Koopmans' Theorem 

What is the physical meaning of the orbital energies e m in Eq. (4.1.28)? This 
was first discussed by the Dutch-American scientist Tjalling C. Koopmans 
(Physica 1934, 1, 1 04.) in his PhD thesis. (After his Ph.D. with Hans Kramers, 
Koopmans began a scientific career in economics. He was awarded the 1975 
Nobel Prize in economics "for his contribution to the optimal allocation of 
resources"). 
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Let us consider the ionization energy of a molecule in the Hartree-Fock 
approximation: 



Note that we calculated the HF energy for two systems, the system N 
electrons and the ionized system of N — 1 electrons. HF approximation will in 
general give different orbitals to the two systems, hence the notation: (p n for the 



that the two sets of orbitals are identical. Of course they are not, but it is 
known that often they are similar so we neglect their difference. Actually, 
when the system is very large and is homogeneous (repeats itself), like an 
infinite crystalline solid, this assumption is expected to be exact because the 
orbitals are spread out on the entire system and therefore removal of just one 
out of an infinite number cannot make a difference. For molecules this 
assumption is a severe approximation. Nevertheless, under this 
approximation we see that all the one body terms cancel except the last and a 
large cancelation of two body terms takes place as well. Only the terms which 
involve the removed orbital stay, these include direct and exchange terms: 



IP, 



HF — 



E H f[4>i> ->4>n-i] - E HF [(p 1 ,...,(p N ] 



(5.2.1) 



N electron system and n for the ionized system. Now Koopmans assumed 



N 




(5.2.2) 



1 N 

+ o / C(mN\mN) - (mN\Nm)) 




Woring out the expression gives: 




(5.2.3) 
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The conclusion: e Ne is the HF approximation to the ionization energy. 
Similarly, £n b -i approximates the next ionization energy etc. 

The flaw in this "theorem" is the neglect of orbital relaxation. In the next 
section we will give a generalized formulation of Koopmans' theorem which 
is exact. In the section after that we discuss the homogeneous electron gas 
which is a system for which orbital relaxation does not exist. 




BNL* BNL* HF HF LSDA LSDA B3LYP B3LYP SAOP KS BNL* BNL* HF HF LSDA LSDA B3LYP B3LYP SAOP KS 

SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. 




BNL* BNL* HF HF LSDA LSDA B3LYP B3LYP SAOP KS BNL* BNL* HF HF LSDA LSDA B3LYP B3LYP SAOP KS 

SCF/TD Kpop. SCF/TD Kpop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. SCF/TD Koop. 



Figure 4: Estimation of ionization energies in HF, approximate and accurate DFT methods. 
Two approaches are given for each method, one based on Koopmans' approach, using the 
orbital energy of the neutral. One can see that Koopmans' HF orbital energies usually 
overestimate the IPs by ~2 eV. Taken from ref. [12]. 

Examples of the performance of Koopmans' theorem within Hartree Fock and 
some DFT brands are give in the Figure 1-1. In the figure the first IP is always 
the first left hand (brown) bar. Also shown are calculations for the IP's using 
the ASCF approach, where the IP is simply the difference between the HF 
energies of the cation and the neutral. In HF theory the Koopmans' approach 
for the first IP are off by ~1 eV ASCF has errors tend to be a bit larger. The 
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usual brands of DFT, such as LSD A, B3LYP have large Koopmans' errors but 
small ASCF errors. The brand called BNL has small errors in both respects. 

C. Fractional occupation numbers , the HF orbital 
functional and the generalized Koopmans ' 
theorem 

We now make a fundamental generalization of the Hartree-Fock energy. 
Instead of viewing it as an expectation value of a determinantal wave 
function, we view it as a new fundamental concept: an orbital junctional. We 
take Eq. (4.1.9) and write it as a functional of "all" orbitals and occupation 
numbers: 

£hf[0i.02.-;A./2.-] 



Here the functional depends on an (in principle) infinite set of orthonormal 
orbitals and their corresponding occupation numbers. Since electrons are 
fermions each occupation number is limited to the unit interval, i.e. < f i < 1 
(since you cannot have negative occupation and you cannot have more than 1 
electron in a given orbital because of the Pauli principle). Let us minimize this 
functional for N electrons. What we mean by this is that we minimize this 
energy E HF with respect to the orbitals n (r) and the occupation numbers f n 
under the constraints that the orbitals are orthonormal and the occupation 
numbers are non-negative, not greater than 1 and that they sum up to N. Thus 
the Lagrangian for this constraint minimization is 
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(5.2.4) 



n 



nm 



W[0i.02<-;/i./2<-] 

= E HF [(p 1 ,(p 2 ,...;f 1 ,f 2 ,...] -2 j / n ^n[<0nl0n) ~ 1] 

n (5.2.5) 

The derivative with respect to the orbitals give the analogous HF equations, 
which are similar to Eqs. (4.1.28) with an important difference: all sums on N 
orbitals are replaced by weighted sums: £n=i^n -* Unfn^n- m essence, the 
usual HF theory is the constant occupation numbers f n = |J n > iV ' ^ 6 
new equations are: 

^ll0n> + V H \<p n ) + K H \(p n ) = £n\(pn) 
1 fm\<Pm(x')\ 2S 



\ m / 
/ S ,\ zVa M f/m0m(*')**K»l) , 

[x K ^ = -e 2 2 j 0mW J i^TTT] dx 



(5.2.6) 



Now, let's discuss the equations obtained by demanding variational behavior 
with respect to f n . Since we have a constraint that f n £ [0,1], we must 
differentiate between three cases: the f n = 0, f n = 1 and fractional (0 < f n < 1) 
cases. To derive with respect to f n is allowed only when you are not at 
constraint boundaries. So in the fractional case an arbitrary infinitesimal 
change in f n is indeed meaningful and the derivative of L with respect to f n is 
equal to zero at the variationally optimal point. This gives: 

X = = yPn\hi\4>n) +2 J /m(< nm l nm > _ (nm\mn)) = e n (5.2.7) 

m 

This equation shows that all fractionally occupated orbitals must be of the 
same energy - X- We can divide all orbitals to be full (f n = 1) or empty 
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(fn = 0) an d partially occupied but these latter orbitals all have the same 
orbital energy. We have seen that for the Hartree-Fock solution to be stable, 
only the lowest orbitals can be occupied. Thus the picture that emerges is: 
there are N — j fully occupied orbitals. Then there can be several degenerate 
orbitals with fractional occupation that adds up to j (assuming N is integer). 
The rest of the orbitals are unoccied.. 

Eq. (5.2.7) is now an exact formulation of Koopmans' theorem. Suppose we 
have a system with slightly less than an integer number of electrons: N — r\ 
where r\ is a small fraction. Then: 

E HF [N] - E HF [N - V ] = e H {N)r] + 0(r) 2 ) (5.2.8) 

Where e H (iV) is the energy of the highest (can be partially) occupied orbital 
for the iV-eletrons system. Thus we see that if we "slightly" ionize the system 
Koopmans' theorem holds exactly (in a sense, the change is so small that there 
will be no orbital relxation). Of course, in real molecules there is no such 
thing as a fractional electron. But still, in terms of the orbital HF theory there is. 

This new concept of an orbital functional has allowed us to consider an exact 
and generalized version of Koopmans' theorem. We will see in later chapters 
that orbital functionals play an important role in advanced approaches to 
DFT. 

D . Hartree-Fock for the homogeneous electron gas 
i. Hartree-Fock orbitals and orbital energies of HEG 

Let us now apply the Hartree-Fock theory to an important system which is a 

model for the valence electrons of a simple metal such as sodium. This is once 

again the homogeneous electron, gas of N electrons in a cubic box of volume 

V under periodic boundary conditions. Note that we are imposing on our 
Electron Density Functional Theory Page 121 

© Roi Baer 



electron gas to be uniform, although it is not necessarily known that this is the 
lowest energy solution. In Hartree-Fock theory for example there are non- 
homogeneous (symmetry broken) states which lead to lower energy. We will 
not study these here however. The density is n = N/V. For an interacting 
system one must have charge neutrality so a static positive charge of the same 
density as the electrons is smeared in the box. This positive charge is called 
"Jellium". The e-e direct term, the Jellium self energy (positive-positive 
interaction) and the electron-Jellium energy must all cancel each other: 

Jee + / vnd 3 r + j pp 

e 2 n 2 rr d 3 rd 3 r' 2 2 ff ^ r ^' r ' 

II "F^T ~ eU JJ Y^¥\ (5.2.9) 



2 

e 2 n 2 ff d 3 rd 3 r' 



So in the Hartree-Fock approach, all that is left are the kinetic energy and the 
exchange energies: 

e HF (n) = t{n) + e x (n) (5.2.10) 

The HF equations are: 

k V 2 cp k (r) + K<p k {r) = e k <p k {r) (5.2.11) 



2m e 



The solutions must be plane waves (since all points are equivalent). Thus: 

h 2 k 2 

+ V (5-2.12) 

And: 

K<p k = X^4> k (5.2.13) 

Let us calculate the exchange eigenvalues X^ 1 . Assume double occupancy of 
all orbitals. Then: 
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K<p k (j) = 



2e 2 
V 3 / 2 



J 



g—iq-r' gik-r 



■ d 3 r' 
\r — r'\ 

e i(k-q)r' 



{q\e q <e F } 
2e 2 r e i\K-q>-r 

> e"? r -d 3 r' 

l—i J \r-r'\ 



y3/2 



e iqr 

{q\e q <e F } 

2e 2 v f e £(fc-<?),(:r '~ r) 



(5.2.14) 



y3/2 



\r' — r\ 



■d 3 r' 



{q\e q <e F } 

Here the summation symbol with {q\e q < £ F ] means that we sum over all 

e iwx 

momentum states q with e q < e F . We now discuss the evaluation of J — — d 3 x. 

\X\ 

Let's assume the box is a sphere of radius R and add a damping factor rj > 
(which can be set to zero after the calculation is done). For a finite volume 



x 



■e vx d 3 x 



= [ e""* [ 
Jo ■'o 



n gVwx cos 8 



-L 



R g{iw-i])x _ g(-iw-r))x 





2?r 
iw 

2n 
iw 



iwx 2 



sin 6 d6 2nx 2 dx 



2nx 2 dx 



e (iw-r))R _ ^ e (-iw-7])R _ ^ 

+ 



(iw — 77) (iw + 77) 



(iw — 77) (iw + 77) 



An 



(1 - e^ R cos(wR)) - —e'^sm 



1 
w 



in(wfl)] 



w 2 + r\ 2 

Take the limit R -> 00 and offer that the limit 77 -> 0. This will lead to the result 



477" 



— . For a finite R, take the 77 -> limit first: 



/ 



iwx 4n 

d s x = — 1 



[1 — cos(wi?) ft finite 
1 i? -> ex. 



The result of inifinte R becomes undefined in the second case. So we use the 
first. Thus we write: 
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2 v -1 4ne 2 1 

{q|e,<e F } 



Where: 



i 2 y 1 4?re2 

^ = _ K Z. (fc - q) 2 (5.2.16) 

{q|e q <£ F } 

Now, to proceed, we want to assume that the momentum states get filled up 
in just the same way as they did in the Thomas-Fermi approach for a non- 
interacting gas. This will be valid if the orbital energies e k are ascending 
functions of k, which according to Eq. (5.2.12) is valid when — is an increasing 

function of k as well. We can only know that however, after we evaluate the 
summation in Eq. (5.2.16)... What we can do, is work in the spirit of the self 
consistent field approach. We will assume that e k increasing functions of k 
then sum Eq. (5.2.16) and check, for self consistency, that — is an increasing 

function of k. Under this plan, we replace in Eq. (5.2.16) the summation over 
all q with e q < e F by a summation over all q with q < k F , which can be 
approximated by the an integral (Eq. (2.1.8)), we find: 

To evaluate this integral we pass over to spherical coordinates, k, 6 and 0: 

1 2 f kF r n sinOdO 2 ^ (5 218) 

%k 71 Jo Jo q 2 + k 2 — IkqcosO^ ^ 

The factor of 2n is a consequence of the integral over the angle (the 
integrand is independent of 0). The integral over 9 is can be performed after 

change of variables x = sin 0, noting that f \ -= — ^ = ( q +k — 2k 1 so: 

& ' & J -lq 2 +k 2 -2kqx \ -2kq ) _ x 
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1 2 r kF k + q 2k F (k F \ 

r k = ~^L ln k^i qdq = -^r F {T) 

We changed variable to x = - and gave the integral over x a new name: 



(5.2.19) 



Fix) 



t rX 



In 



1 +x 



1-x 



xdx = 1 + 



In 



x + 1 



x- 1 



(5.2.20) 



This function is shown in Figure V-5. It is monotonically increasing from zero 
reaching 2 when x -* oo. 



5 1.0 




Figure V-5: The function FQt) defined in Eq. (5.2.20) 
The HF orbital energies are now 

h 2 k 2 e 2 2k, 



2m e 4ne 



2kp /kp\ 
o n V k J 



(5.2.21) 



Clearly, as F(x) is increasing as a function of x, it is decreasing with k, yet it 
comes into the orbital energy expression with a minus sign so overall e k is a 
increasing function of k. A plot of e k , in atomic units is given in Figure V-6. It 
is seen that indeed e k is increasing monotonically with k. 
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Figure V-6: The Hartree-Fock orbital energies e k for various values of k F (all in atomic 
units). For any given k F the function e k is increasing with k. 

Exercise: Determine the Hartree-Fock density matrix of the HEG. 

We have shown that the Hartree-Fock orbitals are plane waves with wave 
vectors less than k F . Thus: 




(5.2.22) 



k<k F 



Denoting r — r' = s we find that p depends only on s and so: 



p(s) = 2 




ik s j3 



d 3 k = 



(2tt) 3 J J. 



2 r k F r 1 



dx 2nk 2 dk 




2 4?r [ sk F . 



1 8?r 



(sin x — x cos x) 




(2?r) 3 s 3 




Since fcp = 3n 2 n we find: 




sk 



(5.2.24) 
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We know that p(r, r) = p(0) = n. A plot of the DM is plotted here: 

p(s)/n 
1.0 Jk 




ii. The density of states of the HEG 

An interesting quantity, determining many of the properties of the substance 
is the density of single-particle states (DOS) D(e F ) at the Fermi level. One can 
calculate and obtain: 

£>0) = j "tide - e k )4nk 2 dk 



V C ep dk 
= 7Tv S ^ e ~ ^nk(e k y—de k (5.2.25) 
{2ny J de k 



V 
2t? 



■6{e F -e)k{eYkXe) 



Here we need to invert the relation e k to k(e). Clearly one can write this as: 

(5.2.26) 



D(e k )=^e(6 F -6 k )k 2 (e'(k)) 



-l 



Thus, to obtain the DOS we need to take the inverse of the derivative e' k . Let 

k 2 j 

us look more closely into that. Note that for free electrons = A /2 / u 3 e thus 

D (e) oc yfe . In particular, the DOS at the Fermi level is finite, a typical situation 
for metals, which is what the HEG is. Now, let us look what happens 
according to HF theory. First, let us look numerically-graphically at e k for 
some arbitrary definite value k F = 1. This is shown in Figure V-7 (left). It 
seems perfectly OK. However, when we plot the derivative e' k , one notices a 
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divergent behavior at k -> k F . Indeed, taking the derivative of Eq. (5.2.21) we 
obtain, using F' (x) = \ (-1 + \ (x + i) In ||±1|): 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

k k 




(5.2.27) 



Indeed, as k -» /c F there is a logarithmic divergence of this expression. Since 
the DOS is proportional to the inverse of e' k , this means that at e F the DOS is 
zero! This prediction by Hartree-Fock theory can easily be checked 
experimentally. For example, electronic conduction of metals is high. Since the 
conduction depends on the availability of electrons at the Fermi level we see 
that HF predicts small metallic conductance, failing miserably... 
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iii. Stability of the HEG in Hartree-Fock theory 

The density matrix assuming unpolarized gas (all electrons are spin-paired) is 
given by Eq. (5.2.24). Now, we can calculate the total exchange energy. From 
Eq. (4.1.17): 

4jj \r-r'\ 

V f v 1/3 /sin k F R - k F R cos k F R\ 2 
"ij. ( '-)**RdR= (5.2.28, 

Vkt f kpv /sinx — x cosx\ 2 



Vkp r Kpv /sinx — xcosxy 
7T 3 J V x 3 / 



xdx 



T , . , rx /sinx-xcosx\ , cos 2x+2x sin 2x-l-2x z 1 , 

Notice that L ; xdx = V - . Thus for x 

J \ x 3 ) 8x 4 4 

oo which is the limit of large volume, we have: 



f /sinx — x cos x\ 1 „ ^ 

K — ? — )^=? (5 - z29) 

Thus for large volume: 

K = --f^ (5.2.30) 

47T 3 



From (2.1.15) n = so: 



47T 

The exchange energy per particle is: 



3kp 



K = - 3 lL N (5.2.31) 



e x = - 1 ±=-C x n^ (5.2.32) 



3 /3X 1 / 3 

With Cjf = - (-) — 0.73856. Let us recall the kinetic energy per particle for 
the HEG (Eq. (2.1.18)): 
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e HF = t + e x = C TF n 2 ' 3 - C x nV* (5.2.33) 

With C TF = 2.871. We see that for high density the energy is primarily kinetic 
and rises with density. For low density the energy drops as the density is 
increased. Thus there is an equilibrium point. This can be seen in Figure V-8. 




0.000 0.005 0.010 0.015 0.020 0.025 



n 

Figure V-8: The total energy, according to the Hartree-Fock approximation, of Jellium, as a 
function of the density n. All quantities - in atomic units. 

We see that Jellium is most stable at density n* calculated from: 

e' HF = \c TP nV* -\c x n~W = -> n* = {^J (5.2.34) 

This gives: n* = 0.002 Iuq 3 and e* = -0.048£ h = -1.29 eV. 

For comparison, sodium, which is a monovalent metal with properties similar 
to Jellium has a total energy per atom (valence electron) of: e Na = — 1.13 eV 
and n Na = 0.0038 . 
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VI . The Hohenberg-Kohn density theory 

In view of the poor predictions of chemical bonds and molecular properties 
afforded by HF approximation and the high numerical price of wave function 
approaches, it is beneficial to seek out methods that circumvent the need to 
represent the many-body electronic wavefunction. We studied in detail two 
theories. One was based on the density but had no real rigorous basis. The 
other was a method that assumed the electronic wave function is of the form 
applicable only for non-interacting electrons. We now want to describe a 
rigorous method that combines ideas of this type in a new way which is both 
rigorous and leads to very accurate approximations. 

A. The first HK theorem 

In electronic structure theory the Hamiltonian is given as: 

H = f + U + j v(r)n(r)d 3 r (6.1.1) 

Where all symbols have been defined in XXX. The identity of the molecular 
system is captured in the external potential v(r). The other terms are 
"universal" i.e. the same for all molecules. In DFT they have a special symbol: 

F = f + U. (6.1.2) 
Different Hamiltonians differ only by their external potential-density term: 

H — H' = j (v(r) - v'(r))n(r)d 3 r (6.1.3) 

This observation has a fundamental implication. Suppose we have two 
electronic wave functions and *¥ 2 of N e electrons which have the same 
density i.e. for all r: 

CFjnCr)!^) = 0F 2 |fi(r)|Y 2 > = n(r). (6.1.4) 

Now, if 
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(¥ 1 |H|¥ 1 )<{¥ 2 |H|¥ 2 ) (6.1.5) 

Then: 

( x P 1 \Pf¥ 1 )<(W 2 \Pf¥ 2 ) (6.1.6) 
Thus the inequality is independent of the position of the nuclei: only the wave 
functions affect it through the universal operator F. The external potential 
term drops out because both wave functions have the same density. 

An interesting natural conclusions is that if Eq. (6.1.5) holds for one 
Hamiltonian the it holds for all Hamiltonians:: 

(¥ 1 |//'|^ 1 )<(¥ 2 |H'|¥ 2 ) (6.1.7) 
This fact will now be used to prove the first theorem of DFT, due to 
Hohenberg and Kohn: 

Theorem (Hohenberg-Kohn): When two Hamiltonians differing only by a 
single particle potential term H — H' = f[v(r) — v'(r)]n(r)d 3 r have non 
degenerate ground states which integrate to the same density then these 
Hamiltonians are identical up to a constant (i.e. v(r) = v'(r) + const). 

Proof of the HK theorem: Assume otherwise: *P is the GS of H and W that of 
H', both wave functions assumed real and have the same expectation values 
for the density at all points in space. The variational principle for dictates 
< (¥'|H| l F') which, as discussed in Eq. (6.1.7) means that < 
holds as well. But this latter inequality contradicts the variational 
principle for H 7 ' as the ground state of H', unless it differs from H by at most a 
constant. ♦ 

This theorem allows one to think of the potential as a functional dependent on 
the density. Thus, in addition to the "usual" 

v H ^x/j ^ n (6.1.8) 
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We now have: 

n -* v -* H -* xp (6.1.9) 

Thus "everything" about the molecule (all its properties) is in the above sense 
a functional of the ground-state density. 

A generalization of the theorem, giving an inequality is: 

Theorem: If rij (r) is the density of the non-degenerate N-particle groundstate 
of H t = F + / Vi(r)n(r)d 3 r, where i = 1,2. Then, denoting AX = X x — X 2 , we 
have: 

An gt => j An(r)Av(r)d 3 r < (6.1.10) 

Proof: Suppose An 0. Then W 2 and because of non-degeneracy the 
following inequality is strict: 

W^iK) < (^2 |^i |^ 2 ) = (^2 1^2 ^2) + J biO) - u 2 (r)]n 2 (r)d 3 r (6.1.11) 
Denoting E t = (^1^1^), we find: 

El <E 2+ jA M r^r (6.U2) 
And exchanging the indices 1^2: 

E 2 <E 1 - j Ai?(r)n 1 (r)d 3 r (6.1.13) 
Adding the two inequalities and cleaning up gives Eq. (6.1.10), QED* 
B. The HK functional 

Since "everything" is a functional of the density, we can assert that the 
ground-state kinetic energy is a functional T[n] of the density and so is the 
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electron repulsion energy U[n]. We can thus define the Hohenberg-Kohn 
functional of v-rep densities as follows: 

F HK [n] = T[n] + U[n] = (W gs [n]\f + 0\W gs [n\) (6.2.1) 

Where ^^[n] is a ground-state wave function with density n. F HK [n] is a 
universal functional, it is not limited to any particular molecular system. It is 
valid for all systems. We have of course no practical way to calculate F HK [n] in 
general. 

Even the domain of definition of F HK is difficult to characterize. In fact, and 
perhaps unexpectedly, this domain is not even convex, as discovered by Levy 
and Perdew. 



Figure VI-1: A convex set is a set of points such that if A and B are in the set then any point 
C on the striaight line joining A and B is in the set as well. The left set is convex while the 
right set is not. 

This means, that if n and n n , 2 are two densities which are in the domain, it is 
not guaranteed that the convex sum n g = (cos 2 9)n + (sin 2 0)n n / 2 (where 6 

is some parameter in the range [o<~]) is in it. Levy and Perdew showed an 

example where convexity fails by considering the case where both are 
densities are densities of two degenerate eigenstates of the same system (we 
discuss this in more detail below). Technically, we say that the domain of 
definition is not convex. However, suppose we do have three densities n , 
n n / 2 and rig, where the latter is the convex sum of the two former densities. 
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Suppose further that all three belong to the domain of definition of F HK . Then 
F HK [n g ] < (cos 2 9)F HK [n ] + (sin 2 0)F HK [n n / 2 ], i.e. F HK is a convex functional 
of the density (we discussed this concept in Chapter XXX). This is seen as a 
result of the variational theorem. Suppose x ¥ g is the ground-state wave 
function, v g the potential and H g = f + U + V g the Hamiltonian 
corresponding to n g then: 

(V |// e |Ve)<(V o |# |%) 



(Ve|^eK><K/2|^eK/2> 

From this we have: 

F HK [n g ] = (*F e |/? e |*F e ) - f n g (r)v g (r)d 3 r 

< (cos 2 6) [(¥ |#eh>) - f n Q {r)v g {r)d 



(6.2.2) 



+ (sin 2 6) 



K/iWbK/i) ~ j n n/2 (,r)v g (,r)d 3 r 



(6.2.3) 



= (cos 2 6)F HK [n ] + (sin 2 0)F HK [n n/2 ] 

We will see that the convextity property of F HK [n] has an important 
implication for the variational property of density functional theory that we 
discuss in the next section. 

C. Minimum principle for density functional theory 

The second HK theorem establishes a minimum principle involving the 
density and it can be used to "find" the density without direct reference to the 
concept of a "wave function". 

Given a potential v(r), consider the following functional for N electrons: 



[n] + j v(r)n 



E v [n] = F HK [n] + I v(r)n(r)d 3 r (6.3.1) 
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The HK theorem II states that the density n v which is the ground-state of v(r) 
minimizes this functional. For n v we know of course that: 

E v [n v ] = (V gs [n v ]\H\V gs [n v ]) = E gs (6.3.2) 

Where H = f + V + U. HK theorem II states that for any other v-rep density n' 
of N electrons: 

E v [n v ] <E v [n'] (6.3.3) 

The proof is an immediate consequence of tha quantum mechanical 
variational principle: 

E v [n'] = (% s [n']\R\V gs [n']) > (V gs [n v ]\H\V gs [n v ]) = E gs (6.3.4) 

This theorem allows one to speak of a "minimum-principle" concerning the 
density. If we have an approximate E v [n] we can find an approximation to the 
ground-state density n v , simply by finding the minimum. 

The fact that F HK [n] is convex (in the limited sense, since there is still the issue 
of domain of definition discussed below), as proved in the previous section, 
implies that E v [n] is convex as well, since by adding a linear term to a 
function one cannot change its convex/non-convex character. The convexity is 
desirable since it assures that the minimum is not only global, but that there 
are no local minima as well, since a convex function (and functional) have no 
truly local minima. There is still the problem of the convexity of the domain and 
this spoils this useful conclusion. Below we show how to construct a 
functional with a convex domain. But then, we lose the convex property of the 
functional. 

We now derive the basic equation of DFT (ignoring for the time being the 
problem of convexity of the domain). We need to minimize E v [n] under the 
constraint that J n(r)d 3 r = N e . Thus we need to minimize the Lagrangian: 
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L v ,N P [n\ = E v [n] - y. 



(r)d 3 r - N e 



(6.3.5) 



Where the value of \i is varied until the constraint is respected. The final 
minimizing density obeys: 

8L ViN [n] 8E v [n] 8F HK [n] 
= - e = - -\i = - + v(r) - n (6.3.6) 
dn(r) on(r) on(r) 

This is the basic equation of DFT It has no direct mention of the wave- 
function. Once we find an approximation for E v [n] we can get an 
approximation for n v (r) from this equation. 

The two theorems of HK put some rigor into the Thomas-Fermi 
approximation. In this KH theory E TF [n] is an approximation to E v [n] and the 
TF equations are an application of (6.3.5). Still, we know that TF theory is very 
poor for chemistry. This means that despite the added rigor, the TF 
approximation is too cumbersome for quantum chemistry. 

D.An interesting observation on the variational 
principle of non-interacting electrons 

Consider a system of N non-interacting particles in a potential v(r). Usually 
we may assume that the ground state of this system is a Slater wave function 
O = det[0 x ... W ]. The variational theorem states that their ground state 
energy is given by minimizing: 

N N 

E[v] = ^(0 n |f! + v\cp n ) - £ e n «0 n |0 n > - 1) (6.4.1) 

71=1 n-1 

Where T 1 = — - V 2 , e n are Lagrange multipliers imposing the unit norm of 
each orbital. Minimizing leads to the equations: 

(f 1 + v)0 n = f n n (6.4.2) 
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After solving for the N lowest energy eigenstates the energy is E[v] = £n=i e n< 
the kinetic energy is T[v] = En=i(0n|7i|0n) the density is 
n(r) = 2Lil0 n (r)l 2 . 

Now consider a different problem. Suppose the density is given and one is 
required to find the system with this density that has the minimum kinetic 
energy T. This leads to the Lagrangian: 

n In ' 

T[n] = ^(0 n |f 1 |0 n ) + j v(r) ^J0 n (r)| 2 - n(r) 

(6.4.3) 

N 

-^^«0nl0n)-l) 

n=l 

Now v(r) are Lagrange multipliers and must be searched for in order that the 
constraint n(r) = En=il0n( r )l 2 be fulfilled. After minimizing we obtain the 
equations: 

(f 1 + v)0 n = e n n (6.4.4) 

This equation rises from the attempt to compute T[n] and is the same as Eq. 
(6.4.2) which rose as when v was given and E[v] was calculated. This shows 
that minimizing the kinetic energy under a given density invokes "the 
same"equations as minimizing the energy when the potential is given. This 
fact will be used in the Kohn-Sham method. 

E. The set of V-representabile densities 

i. V-rep densities correspond to ground states wave function of some 
potential well 

The HK theorem shows that the ground state density of a system uniquely 
determines the one body potential. This is a uniqueness statement: there is at 
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most one potential associated with a density. An interesting twist is the 
reverse question: what are the conditions that a given density is the ground 
state density (GSD) of some system? Of course there are some preliminary 
conditions on the density: it must be non-negative and it must integrate to a 
positive integer: 



But in general, we have no good criterion for deciding whether a given 
density n(r) is a GSD of some potential v(r). Densities which are GSDs of a 
potential are called "v-representable". In short v-rep. 

When we say that everything is a functional of the density we mean 
everything is a functional of a v-rep density. 

ii. Some non-v-representabilit issues 

We have seen in XXX that ground state wave functions of single particles is 
nodeless. A corollary from the above analysis is that a density of one particle 
with a node is not v-representable. 

However, the density does not need to actually develop a zero for the density 
to be non-vrep. Consider the example by Englisch and Englisch: 



n(r) > 




(6.5.1) 



n(x) oc (1 + (x 2 )") 2 e- 



1 1 




If ip{x) is a wavefunction then the potential is given by v(x) 



lxp"(x) 
2 1p(x) 



Using 



the above form for a = - , 



we find that the potential is infinite at the origin: 
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Figure VI-2: Examples, following Englisch and Englisch, of a non-vrep density (in purple) 
and the corresponding potential which is singular. The left panel is a ID example which 

the right panel is the radial 3D example. 

In 3 dimensions simply replace x by r and i//'(x) by V 2 i/;(r) = ^(r 2 xjj(r) s ) 

! (r 2 Tp(r))" ! 

thus: v(r) = - . : -. The first term of the potential is seen to 

v J 2 r 2 xp(r) r 2 r 

dominated by a term going as — when r -> 0. This is a centrifugal barrier with 

■6 = 1. Subtracting this barrier exposes the bare potential, shown in Figure 
VI-2. 



iii. The set of v-rep densities of a given electron number is not convex 

Degenerate Hamiltonians can generate non-vrep densities quiet easily. Thus, 
non-vrep densities are much more abundant than one may suspect. Let us see 
why there is a problem, in the following analysis due to Levy and Lieb 
(developed seperately at the same time more or less). 

Suppose H = f + / v(r)n(r)d 3 r + U is a Hamiltonian with degenerate ground 
states, of energy E and full degeneracy Q. Thus, H^, i = 1, ...,Q are Q ground 
state wave functions with (Vj)^-) = 5^-. There are infinitely many ways to 
define the We select one arbitrarily. Denote: n £ (r) = (W t \n(r) \Wi). Then, 
consider the density built as a convex sum of these degenerate-state densities: 
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Q Q 
n(r) = ^ qnjCr) c £ > ^ c £ = 1 (6.5.2) 
t=i t=i 

We now show that this density it is usually non-v-representable. Of course, it 

may happen (in rare cases) that a linear combination of the degenerate wave 

functions ¥ = £? =1 yields n(r) = 0F|n(r)| l I'>. This case however is not 

ordinary and we consider the cases where this does not happen. We thus 

proceed to show that n(r) cannot be the groundstate density of any other 

Hamiltonian as well. We do this by contradiction and assume existence of a 

wavefunction W which is the GS of some H' and such that ( l F'|n(r)| l l / ') = 

n(r). The variational principle states > for each i. 

Multiply by the positive c £ (so the inequality is not spoiled) and sum over i, 

using £?c £ = i, and obtain > Now use the same 

reasoning that led to Eq. (6.1.7) and replace H' by H, obtaining: 

Ci^H^i) > (V\h\V). However, (W^H^) = E and so this leads to 

< E which condtradicts of the variational principle for Wi. Hence, 

H*' cannot be the ground-state of any Hamiltoian. 



Example: The above theorem is general and holds for any system 
of particles. We consider non-interacting electrons. Consider the 
density of non-interacting electrons in the potential well created 
by a a Lithium nucleus. The ground state is 4-fold degenerate 
*Fi = |lsls2s|, ¥ 2 = \lsls2p x \, V 3 = |lsls2p y |, V 4 = \lsls2p z \ 

r r r 

, e~ r e~2( r -2) e~2~x e~2y 

Wlth ^ = _ ^ = __ ^ = _ / ^ = __ ^ 2pz = 

_r 

j==> The 4 densities are: n x = 2n ls + n 2s , n 2 = 2n ls + n 2Px , 
n 3 = 2n ls + n 2pyl n 4 = 2n ls + n 2pz . The average is given by: 

lV e- 2r e~ r 

n(r) = 4 Z, n£(r) = 2 + 128^ ((r " 2)2 + 

t=i 
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This density is spherically symmetric and is plotted below, 
together with the density n x (r) for comparison: 




Figure VI-3: Left panel: The average density of n(r) plotted together with 
the almost identical density n x (r) orbital. Both densities seem almost 
identical yet the latter is v-rep while the former not. Right Panel: The 
difference between the two densities 



This density is very naive and there is no visible indications 
suggesting this is not (non-interacting) v-representable. 

We ave found that such convex sum densities are not v-representable in the 
sense that there is no Hamiltonian for which a ground-state yields the density. 
However, clearly these densities are associated with the potential v(r) which 
created the degenerate states. One can extend the Hohenberg-Kohn functional 
definition to deal with these convex sum densities in the following way. Given 
a density which is not v-representable in the usual way we now call it non- 
pure-state v-representable. Then we will say that such a density is ensemble- 
v-representable and assume it is associated with the potential v(r) from which 
the degenerate ground-states ^ were formed. Then for such a density the HK 
functional is written as: 

Q Q 
F HK [n] = £ c t {%\f + Df¥ t ) n(r) = £ Q<^|^(r)|^> (6.5.3) 
t=i t=i 

Here, we implicitly made the identification n -> v -> q. The HK uniqueness 
theorem guaranteeing that only one potential v can form a given n can be 
extended also to this case (i.e. that n -> v is meaningful) still holds. 
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Furthermore, the minimum principle, i.e. that E v [n] with the more general F HK 
obtains its minimum at the density corresponding to v(r) remains true. 



F . Levy-Lieb generalization of the HK functional 

The minimum principle of HK is of crucial importance for density functional 
theory. Yet, the basic equation (Eq. (6.3.6)) which is the basis for development 
of DFT into a practical approach is somewhat problematic from the 
mathematical point of view. In order to define the functional derivative, one 
must ensure that there is an "open" neighborhood of the ground state density 
n v (r) for which F HK [n] can be defined for any density. As we saw, there is a 
fundamental problem is that of v-representability. The ensemble v- 
representability solved one kind of problems but it is not clear if there may be 
additional classes of non-v-representability. It is not possible to assume that 
any density that is positive and integrates to an integer is a density of a 
ground state wave function of a Hamiltonian. Furthermore, we cannot even 
assume that around a v-representable density there is an "open 
neighborhood" of v-representabile densities. 

In general one cannot waive the possibility that there may be densities which 
are not the ground state densities of any system and still are arbitrarily close 
to any n v (r) of interest. In such a case the functional derivative of F HK is 
formally undefined. 

A way around these problems, developed separately by Levy and by Lieb is 
to formulate a functional F LL [n] defined for any density on the one hand and 
equal to F HK [n] for v-representable densities on the other which still allows 
for the same type of minimum principle as does the HK functional. We first 
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note an important property of F HK [n]. We say that the wave function 
realizes the density n(r) (and symbolize this as *F -» n ) if : 

(¥|n(r)|¥) = n(r) <=> ¥ -> n ( 6 - 5 - 4 ) 
We will now show that for any wavefunction *F which realizes n(r) 

ip^n => F HJf [n] < (V|f + f7| l F) (6.5.5) 

This is a direct result of the variational principle: F HK [n] = E gs [n] — 
f v(r)n(r)d 3 r < (^H^) - f v(r)n(r)d 3 r = -9\*¥) = (v|f + t/|v). 
From this development we see also that the wave-function which realizes this 
minimum is the ground-state wave function and is denoted Iggfri]. One has 
therefore: 

F HK [n]=min{W\f + 0\v) (6.5.6) 

Levy and Lieb decided to use this relation, valid only for v-rep densities as the 
definition of their functional, valid for any density: 

F LL [n] =mm(w\f + D\w) (6.5.7) 



Exercise: Show that for any density there is at least one wave 
function which realizes it 

Solution: For a one-electron case this is trivial, since the density is 
non-negative we can take the wave function as *P(r) = ^n(r). 
For N electrons, just "slice" n(r) to N non-overlapping parts 
each: positive, integrating to 1: (1) n(r) = Efc=i n /e( r ) > (2) n k (r) > 
= 0, (3) 4>i(r)4>j(r) = Sijnj(r) where k (r) = y/n k (r), and (4) 
j n k (r) = 1. This can be done in countless ways (with a good 
sharp knife). Then the determinant = |0i0 2 ■■■ 4>N e \> is a 
wavefunction which creates n(r). 
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The above exercise shows, that the definition in Eq. (6.5.7) makes sense for 
practically all densities. One can now use F LL [n] instead of F HK [n] in Eq. (6.3.1) 
and convert the Hohenberg-Kohn minimum principle to a variational- 
minimum principle. The search for the constrained minimum can be done 
using a Lagrange multiplier approach. We formulate the Lagrangian: 

PM».v].W + m + l«rWWr)m-«<rM*r (6.5.8) 

Given n(r), one now find suitable Lagrange multipliers v(r) such that when 
minimizing F LL with respect to an antisymmetric wave function leads to the 
constraint (^InCr) IV) = n(r) at any r. The search for a minimizing wave 
function will cease when the gradient is zero, so a necessary condition is: 

■j^PuW.nMW^u* = = (f + 0+( v LL (r)n(r)d 3 r)v LL 

°V v=v LL \ J I (6.5.9) 

= Hll^ll 

Thus, when n(r) is given, the Lagrange multiplier function v LL (r) that 
imposes the condition defines a Hamiltonian H LL = f + / v LL (r)n(r)d 3 r + U. 
Note that v LL (r) is determined by the procedure up to a constant. If n(r) is v- 
rep then necessarily, by Kohn's theorem, F LL [n] — F HK [n] and is the 
groundstate wave function. However, we are not in general assured that 
is always the ground state, since n(r) is not necessarily v-rep. Indeed, from 
(6.5.9), all we can say is that H^l must be some eigenstate of H LL . 

Once the minimum is attained, one can take the functional derivative of F LL by 
n(r). As with any constrained minimization, when the derivative of the 
optimized solution is taken with respect to the constraints, one obtains the 
Lagrange multiplied (see (1.11.23)): 
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P LL \V,n,v] = -v LL (r) (6.5.10) 



Sn(r) 

Now, once the constrained minimum F has been found, one can plug F into 
the HK scheme and consider the Lagrangian corresponding to E v = F LL + 
j v(r)n(r)d 3 r: 

L v [n] = F LL + j v(r)n(r)d 3 r — fi J n(r)d 3 r — N (6.5.11) 

If now one searches for the density that minimes this functional one 
immediately obtains: 

= = -%(r) + v(r) - n (6.5.12) 

Showing that the minimum of E v [n] is obtained when the density found 
admits, as Lagrange multiplier v LL (r), the same potential as was give to start 
with v(r), up to a constant. 

Even with the LL functional one problem is still not solved: convexity. While 
the domain of definition is now finally convex, the functional has lost 
convexity. A counterexample for convexity can be given, using the degenerate 
system used above to prove that the domain of definition of the HK functional 
is non-convex. Indeed, for that system, which had a potential v (r) and 
degenerate eigensstates W t (i = 1, ... , Q) all of energy E: 
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LL 



= (V LL \f + U\V LL ) 

i 

> ^ c { E — j c i n i (r)v(r)d 3 r (6.5.13) 

i 

= ^ c t [(Vi|T + + V^j) - j ni(r)v(r)d 3 

i 

= 2^c i F LL [n i ] 



Showing that for these densities the functional is concave (while for the case 
where rij and £t c £ n £ are all v-rep the functional is convex). Therefore overall 
the functional is not convex. 



Once again, it is possible to circumvent this problem by extending the 
definition to ensemble densities. 



Exercise: Discuss the Levy approach for a given density n(r) of a 
system of non-interacting electrons. 

In this case the Hamiltonian is H = t + / v(r)n(r)d 3 r i.e. a 1- 
body Hamiltonian. Suppose ip m/ m = 1,2, ... are the 1-particle 
states which are eigenstates of fi x = 7\ + / v(r)n(r)d 3 r: 

h 2 

= f m 0m( r ) 

Let us suppose orbital energy ordering, so that e x < e 2 < e 3 ... If 
the density is "non-interacting v-rep" then it must be due to the 
Slater wave function D gs = |0i"-0w e | of the lowest energy N e 
orbitals. In this case, the kinetic energy is given by: 

N e 

T LL [n] = £(0 m |f|0 m ) (6-5.15) 

771=1 

And the functional derivative is: 
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-S^r) T iL[n\ = -v(r) (6.5.16) 

If we are given a non-interacting non-v-rep density then we 
might find not a ground state but an excited determinant, for 
example: D ex = |0i •■• 0w e -i0w e +i|- This leads to "holes" in the 
non-interacting system. 



G. The dilation inequality for the HK functional 

We have already seen the the HF functional is convex. More exact properties 
are desirable, so that when one derives approximations to this functional, they 
can perhaps be forced to obey the known exact relations. One approach is to 
use dilation considerations. This will lead to an interesting inequlaity. 

We will use the following symbol convention: 

*P y (r lf r 2 ...,r Ne ) = Y 3Ne/2 V(yr lf yr 2 -,yr Ne ) 

(6.6.1) 

n y (r) = y 3 n(yr) 
It is quite straightforward to show the following relations: 

Exercise 

Derive the following relation 

(¥ y |n(r)|¥ y ) = y 3 ^\n{yr)YV) 

(V Y \f\V Y ) = Y 2 (v\f\v) (6-6-2) 

From the first relation in this equation, one can deduce that if -> n(r) (i.e. T* 
is a many-electron wave function exhibiting the spatial density n(r)) then: 

(^InCr)!^) = y 3 0F|n(>r)| l I'> = y 3 n(yr) = n y (r) ( 6 -6.3) 

One can now also use this relation to get a dilation inequality. Indeed, since 
Vjnjy realizes the density n y (r) (see the first of Eqs. (6.6.2)) we have 
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F HK [n Y ] < (V[n] Y |f + ^|V[n] Y ) = y 2 (V[n]|f \W[n]) + Y^ln^O^in]), thus we 
find: 

FhkM ^ Y 2 T[n] + yU[n] (6.6.4) 

VII. The Kohn-Sham method 

Kohn and Sham noticed that the HK theory is valid for both interacting and 
non-interacting electrons. Now, they ask, what happens if for any system of 
interacting electrons, with density n there is a non-interacting system of the 
same density? It is clear that if such two systems exist they are unique. The 
non-interacting system has one advantage over the interacting system: we can 
find its ground-state rather easily, since the many-body wavefunction is a 
Slater wave function. So, the problem is: how to perform such a mapping. 

A. Non-interacting electrons 

If non-interacting electrons are tractable, let then study their density 
functionals. First, if we are given n(r) we assume it is non-interacting v- 
representable. That is we assume there exists a potentials v s (r) such that the 
ground-state of the resulting non-interacting Hamiltonian H s = T + V s admits 
a ground-state having n(r) as density. We can later use the Levy-Lieb 
approach to patch up the v-representability requirement. 

Let us denote the ground-state wave function of non-interacting electrons that 
have the density n(r) by 5S [n]. The Hohenberg-Kohn functional for non- 
interacting electrons is reduced to just the kinetic energy, i.e. we define: 

T s [n] = (% s [n]\T\% s [n]) (7.1.1) 
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An interesting question is - can we devise a method to compute T s [n]. For 
most densities, the non-interacting wave function is simply a Slater wave 
function. The minimum principle for T s [n] (derived from Eq. (6.5.5)) is: 

ip^n =^ T s [n] < (W\f \W) (7.1.2) 

One should not treat this equation lightly and realize its non-trivial 
consequences: that the ground-state wave function of non-interacting particles 
associated with n(r) minimizes the kinetic energy! Let us see some 
consequence this minimum principle. 

Let 5S [n] be the ground state wave function of the system of non interacting 
electrons realizing the density n(r). From the first equation in Eq. (6.6.2) one 
sees that 5S [n] y realizes n y (r), thus one can plug it into the right hand side 
of Eq. (7.1.2), with n y (r) plugged into T s : 

TsM = KsKP Ksk]> < (%sin] Y \f\% s [n] Y ) (7.1.3) 

Since now the right hand side is an expression of a scaled wave-function, one 
can use Eq. (6.6.2) and obtain: 

T sM < (%s[n] Y \f\% s [n] Y ) = r 2 (O flS [n]|f \% s [n]) = Y 2 T s [n] (7.1.4) 

And so T s [riy] < 7 2 r s [n]. But we can also use this equation to scale n y back to 
n by dilating by 1/y, since: { n y)y Y = n - Then, the same rule applies and we 

get: T s [n] < Y~ 2 Ts[n Y ], and so T s [n y ] > 7 2 r s [n]. We obtained two contradicting 
equations which can agree only if both are reduced to equality. Thus: 

T s [n Y ]= Y 2 T s [n]. (7.1.5) 

This should be compared with the analogous result of Eq. (6.6.4) which is an 
inequality. The dilation effects for non-interacting electrons is obviously much 
simpler. One important corollary of (7.1.5) is that the non-interacting wave 
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functions scale with the density. This can be seen if one remembers that 
5S [n] minimizes the kinetic energy. Therefore any function 0[n] for which 
(0[n]|f |0[n]) equals T s [n] must be the ground state of the non-interacting 
system with density n. We saw that 5S [n] y gives the correct kinetic energy of 
for the system with density n y , Ts[n y ], and therefore, necessarily: 

O flS [n] y = ^sDv] ■ (non — interacting particles) (7.1.6) 

Exercise 

1) We can define a functional called the Hartree energy which is the 
classical electrostatic energy associated with the charge distribution 
n(r): 

- jR^*T*V (7.1.7) 

Prove the following dilation relation: 

E H [n Y ]= Y E H [n] (7.1.8) 

2) Now define the exchange energy functional (see also Eq. (4.1.17) for a 
definition based on the orbitals): 

K[n] = (* gs [n]\0\*g S [n\) - E H [n] (7.1.9) 

Use (7.1.6). Prove: 

K[n Y ] = yK[n] (7.1.10) 

What is the relation between the potential v s (r) for which n(r) is a non- 
interacting ground state density and Vgyir) for which n y (r) is a non- 
interacting ground state? We can use the basic DFTequation to answer this. 

From the basic definition of functional derivation: 
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ST, 



s 



8n(r') 



- lim fefrfrfrj + " r,) l " ^K^)]) (7in) 



The 3D delta-function has the density dilation structure: S(r — r') = 
A 3 5(lr - Ar'), so: 

,. ( T s [y 3 ( n Gr) +J7ff(yr-yr'))] - r 5 [y 3 n(yr)]) , 711 ^ 

M-y - v S y(r') = hm-^ — (7.1.12) 

Then using Eq.: 

f „ ,. K 2 <T s [n(r) + ?/5(r - yr')] - r s [n(r)]) 

v Y ^° V (7.1.13) 
= y 2 - Vs(y')) 

By dilating in the reverse direction we can easily see that: 



v s ,y(r') = Y 2 v s (yr') 



(7.1.14) 



We could have obtained this result directly from the Schrodinger equation. 
Suppose ^({r}) is a many-body eigenfunction of Hamiltonian H = f + K({r}), 
i.e T = (E — vyv. Define a scaled wavefunction: 



3N 



^(ri r «) = y ~^(yri r^w) ( 7 - L15 ) 

Then clearly: 

^({r}) = Y 2 Y--{t^){{yr}) =y 2 (E- V({yr}))*F Y ({r}) (7.1.16) 
And so obeys the S. E. tW y = (E Y — V y y¥y with energy and potential given by: 

V Y ({r}) = Y 2 V{{Yr}) 

(7.1.17) 

E Y = Y 2 E 
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For a ground state of non-interacting the first equation means that v sy (r) — 
Y 2v s(.Y r ) is the one-particle potential for the scaled determinant, and thus for 
the scaled density (since the scaled determinant realizes the scaled density). 

B. Orbitals for the non-interacting electrons 

Consider the ground state wave function for the non-interacting electrons 
5S [n]. Suppose the non-interacting electrons reside in the potential v s (the 
subscript S is in honor of Slater) then (f + %)<t> gs [n] = E s O gs [n]. In most 
cases this wave function is a normalized Slater wave function. The only 
exception might be the existence of a degeneracy and then there may be 
several independent Slater wave functions and O fiS is a linear combination. 
However, one can always assume that he is looking for one of the 
determinants and not for a linear combination (this means that one must be 
careful not to impose possible symmetries of the Hamiltonian). Once one 
does that, we can introduce N e orthormal spin-orbitals (p q (x) q = 1 ...N e , from 
which (t> gs [n] is built. These orbitals are excited states of the potential well in 
which the non-interacting electrons reside. Thus the orbitals must each obey 
the single-electron Schrodinger equation: 

-^ 2 ^W + v s (r)0,OO = e q 4> q (x) (7.1.18) 

The orbitals correspond to the N e lowest eigenvalues e q . The fact that <t> gs 
realizes the density n(r) is expressed as : 

N e 

<r) = ^j0,(*)| 2 (7-1.19) 

q = l 

The non-interacting kinetic energy T s [n] is then: 
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T s [n] = ^ | 0,(x) (-^F 2 W*)dx (7.1.20) 

q=l ^ ' 

When one wants to find the functional derivative of T s with respect to n(r) 
one can turn to basic equation, Eq. (6.3.6) which in case of non-interacting 
electrons becomes: 

■Wrj + vs(r)=n (7.1.21) 

We will see that this equation is important for the method known as the 
Kohn-Sham method. 

Note that the discussion of dilation in the previous subsection can be carried 
on to the orbitals. The only additional information to convey is that the 
orbitals scale as the density and each of the orbital energies scale as the total 
energy: 

e qiY =y 2 e q (7.1.22) 

C. The correlation energy functional: definition 
and some formal properties 

The ground state energy of the an system of electrons in density n(r) can be 
written in terms of the non-interacting (Slater wave function) wave function 

E[n] = {* gs [n]\H\4> gs [n]) + E c [n] (7.2.1) 

This equation is actually a definition of a new density functional, the 
correlation energy functional £" c [n]. If we suppose for the time being that 
E c [n] is known, this expression can be used to define a working procedure 
for DFT known as the Kohn-Sham method. From our studies of the Hartree- 
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Fock theory, we know that the expectation value of H within a determinant 
can be written as: 

(%s[n]\H\% s [n]) = E H [n] + K[n] = U s [n] (7.2.2) 

Where E H is the Hartree energy: 

e 2 ff n(r)n(r') 
E H [n] = y JJ |r _ r , | d 3 rd 3 r' (7.2.3) 

And K is given in terms of the orbitals from which 5S [n] is composed (Eq. 
(4.1.17)): 

K[n] = --{[ d 3 xd 3 x' (7.2.4) 

2 JJ \r-r'\ 

Based on this, we rewrite F HK [n] as: 

PhkM = (V> gs [n]\P\V gs [n]) = (% s [n]\F\<t> gs [n]) + E c [n] (7.2.5) 

This allows us to write Eq. (7.2.5) compactly as: 

F HK [n] = T s [n] + U s [n] + E c [n] (7.2.6) 

Clearly we have also the equivalent equation: 

E c [n] = T[n] - T s [n] + U[n] - U s [n] (7.2.7) 

Physical intuition concerning molecules and solids tells us that E c [n] is a 
small quantity when compared to T s [n] or U s [n]. Thus, it is reasonable to look 
for approximations to this quantity. Approximation to the correlation energy 
functional is the most important issue in DFT It is an open question, still 
being worked upon. 

We shall deal with approximations later. Meanwhile, let us ask what can be 
safely said about the correlation functional. We prove here several important 
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(7.2.8) 



inequalities. First, consider the difference T c [n] = T[n] — T s [n\, the correlation 
kinetic energy. 

Exercise 

Show T c [n] is a positive quantity. This is actually intuitively clear: the 
interacting electrons must have much more complicated "paths" in the 
interacting case because they want to avoid "bumping into" other electrons. 
Anything with more swirls must have higher kinetic energy. 

Solution 

Using the variational theorem: 

T s [n] + V s [n] = (% s [n]\f + V s \<S> gs [n]) < (% s [n]\t + 9 s \W gs [n]) 
= T[n] + V s [n] 

Where W gs [n] is the interacting ground-state wave function determined by n. 
Comparing the two sides we have: 

T s [n] < T[n] (7.2.9) 

Next, we can show that the exchange correlation energy is always negative. 
This comes about from our experience with expectation values of 
determinants: 

F HK [n] + V[n] = E gs < (<P gs [n]\n\t> gs [n]) = T s [n] + V[n] + U s [n] (7.2.10) 
Using Eq. (7.2.6) we find: 

E c [n]<0 (7.2.11) 

Furthermore U c [n] = U[n] — U s [n] is negative as can be seen from the fact that 
E c is negative and T c is positive: 

U c [n] =E c [n]-T c [n] <0 (7.2.12) 

An additional property is the dilating relations. We have proved that 
FhkM - Y 2 T[n] + yU[n] and T s [n Y ] = Y 2 T s [n], U s [n Y ] = Y^sM- Thus: 
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E c[n Y ] = T c [n Y ] + U c [n Y ] < y 2 T c [n\ + yU c [n] (7.2.13) 

One way to proceed is to substitute U c with E c — T c . The other is to substitute 
T c with E c — U c . We obtain from each possibility: 

E c [n Y ] < y(y ~ l)T c [n] + yE c [n] 

(7.2.14) 

E c [n Y ] < Y 2 E c [n] + y(l - y)U c [n] 

And, since T c is always positive and U c negative, we find that: 

E c [n Y ] < yE c [n] y<l 

(7.2.15) 

E c [n Y ] < Y 2 E c [n] y < 1 

Obviously, the second relation is contained in the first (since E c < 0) and so 
only the first relation is important; it can be used to derive complementary 
inequalities. Indeed, applying them for y -> y~ x we find: 

EcWiy\<y- x E c \rx\ y>l (7.2.16) 

This holds for any n, so we stick in (7.2.16) n y instead of n and using the fact 
that (n v ) , = n to obtain: 

E c [n Y ]>yE c [n] y>l (7.2.17) 
D. The Kohn Sham equations 

i. The Kohn-Sham equation from a system of non-interacting problem 

Let us now turn to the issue of determining v s (r), required for the mapping 
between the interacting non-interacting systems. Let us start form the basic 
equation Eq. (6.3.6) which becomes, using Eq. (7.2.6): 

6T s [n] 



Sn(r) 

Where: 



+ v HX (r) + v c {r) + v(r) = /i (7.3.1) 



8U s [n] 

v H x(r) = = v H (r) + v x (r) (7.3.2) 
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is the Hartree-Exchange potential. Notice that: 



v H (r) = 



SE H [n] 
8n(r) 



n(r') 
|r — r'| 



(7.3.3) 



and: 



v x (r) 



SK[n] 
Sn(r) 



(7.3.4) 



Thus, from Eq. (7.1.21) and Eq. (7.3.1) (up to a constant): 



Vs (r) = v (f)+ v HX (r) + v c (r) 



(7.3.5) 



This equation gives us the potential of the non-interacting system. Thus, we 
have made a definite connection between the interacting and non-interacting 
systems. 

Now, an important observation allows us to set up a simple method for 
obtaining the ground-state of an interacting system of electrons. We need to 
find a density that obeys two conditions: 

a) It is the density of the non-interacting electrons so it is the sum of square 
orbitals (Eq. (7.1.19)) that obey the Schrodinger equations (7.1.18) with 
potential v s (r). 

b) The potential v s (r) must be related to the interacting system by Eq. (7.3.5) 

This leads to a simple SCF procedure, called "The Kohn-Sham method" 
reminiscent of the Hartree-Fock algorithm: 

1. Guess n(r). 

2. Build v HX (r), v c (r) from n(r) (Eq. (7.3.3) and (7.3.4)). 

3. Obtain the orbitals (p q (x) from the lowest N e eigenvalues of Eq. (7.1.18). 

4. Compute the density from Eq. (7.1.19) 

5. Redo from step 2 using the new density until you converge - i.e. until 
the density changes no more. 
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When the process convergence, we have the exact ground-state density n v (r). 
It can be used to compute the ground-state energy by plugging it into the 
energy functional of Eq. (6.3.1): 

E v [n\ = T s [n] + j v(r)n(r)d 3 r + U s [n] + E c [n] (7.3.6) 
ii. Systems with partially occupied orbitals 

Partially occupied orbitals arise in the realm of ensemble-DFT, i.e. when the 
density searched for is not "non-interacting v-representable". We still assume 
the density is a sum of orbital densities as follows: 

n(r) =^J<^(*)| 2 , (7 3 7) 

q 

but unlike Eq. (7.1.19) we do not impose exactly N e orbitals. The orbitals are 
almost completely unconstrained, except for the Fermionic condition, that 
their norm must be not be great than 1: 

(4> q \4> q ) = n q < 1, (7.3.8) 
but the total number of electrons is still N e : 

j n(r)d 3 r = ^n q = N e ( 7 3 9 ) 

And the kinetic energy is written similarly to Eq. (7.1.20): 

T s [n] = ^ | 0,(x) (-^-V 2 \<p q {x)dx (7.3.10) 

q e 

We generalized Eqs. (7.1.19) and (7.1.20) any number of orbitals leading to the 
following orbital Lagrangian: 
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L rb[{(pq}. {nq}] 



= E V , 



(7.3.11) 



e q (q = 1,2, ...) and [i are Lagrange multipliers assuring orbital normality and 
total number of particles. We minimize the Lagrangian with respect to the 
orbitals (p q (r) and the occupation numbers (constraint to be non-negative and 
not greater than 1). The constraint minimum procedure thus requires that: 



= ana — = 



8<p q {x) 



Sn 



(7.3.12) 



The first equation leads to the KS equations: 
h 2 

■V 2 (pJx) + v s (r)(pJx) = 6 q (pJx) 



2m, 



(7.3.13) 



For convenience, we order the indexing so that the series of orbital energies is 
ascending: e 1 < e 2 < ••• The orbitals are now eigenfunctions of a Hermitian 
Hamiltonian, and so we can assume they are orthogonal: 

{<t> q \<t> q ,) = n q 8 qq ,, (7.3.14) 

We define as a short-hand notation v s (r) = v(r)+ v HX (r) + v c (r) (see Eq. 
(7.3.5)). By multiplying Eq. (7.3.13) by 4> q (x) integrating on x and summing on 
q we find: 



(7.3.15) 



The second condition for minimum of L orb holds only in cases that the 
minimum is attained with non-integer electron number, n q < 1. The cases 
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n q = 1 or n q = are the boundary of the constraints and the derivative need 
not be zero there. From this second part of the equation, we have: 

e q = ii, <n q <1 (7.3.16) 

This shows that all incompletely occupied orbitals have the same orbital 
energy equal to the chemical potential n, the Lagrange multiplier imposing 
the number of particle constraint of the Hohenberg-Kohn theory (Eq. (6.3.6)). 

Typically all orbitals with e q < \i are fully occupied (n q = 1) and those with 
e q > ii are fully unoccupied (n q = 0). The interacting electron energy is 
obtained using Eqs. (7.3.6) and (7.3.15) as 

B. = 2 e q n q + E HxM + *[»] - J nfrXWrt + r c W)dV (7 . 3 ,7) 
q 

We should note that the development here assumed non-interacting v- 
representability. In cases where this is not valid other occupation rules may 
apply. 

iii. Is the ground state wave function of non-interacting particles always a 
Slater wave function 

The standard KS approach to the non-interacting kinetic energy is by defining 
the functional T s as a minimum principle on the manifold of single 
determinants: 

T s [n] = min (<&|f |o) (7.3.18) 

<J>->n(r) 1 1 v ' 

Setting up a Lagrangian and searching for the constrained minimum yields N 
occupied orbitals. If the density is not v-representable one or more of the low 
energy orbitals may have occupation numbers. 
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Another way to define a non-interacting kinetic energy functional is by an 
extended minimum search over more general wave functions: 

T NI [n] = min Mf \W) (7.3.19) 

Usually this search ends with being a single Slater wave function and for 
these case T NI and T s given are the same. Yet, this may not always be the case. 
Let us assume that we are searching through all W — cos 9 <S> A + sin 9 O b 
where 4 or O s are Slater wave functions. In general, the kinetic energy is 
then: 

T NI [n] = min (0 A |f |O a ) cos 2 9 + sin 29 (0 A |f |O b ) 

^ n(r) (7.3.20) 
+ (d> B |f |4> B )sin 2 9 

Now, if <t> A and O b differ by only one orbital (say (p A = 0# for k = 1, ... , N — 1 ) 
then W is actually a Slater wave function: its orbitals are the N - 1 0^'s and 
then a new orbital 0jj ew = cos 9 (p% + sin 9 0^ is added. Next, if <& A and O b 
differ by two (or more) orbitals the cross term in (7.3.20) is zero and the 
kinetic energy is simply the sum of orbitals kinetic energies with occupation 
numbers given by cos 2 9 and sin 2 9. The orbitals shared by the two 
determinants will have unity occupation number (since cos 2 9 + sin 2 9 = 1) 
while orbitals in the A determinant but not in the B determinant will have 
occupation number cos 2 9 and those in the B determinant but not in the A 
determinant will have occupation number sin 2 9. The orbitals coming out 
from the minimization will all solve a Schrodinger equation with the same 
potential. The 4 odd orbitals (2 from each determinant) will all have the same 
orbital energy equal to pL. When there are more than 2 Slater wavef unctions a 
similar treatment will result and even more orbitals will be degenerate at the 
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chemical potential. Even if n is v-representable this type of wave function can 
arise. 

A third way to define the kinetic energy is as an ensemble average. Instead of 
a wave function, one uses a mixed density matrix D = Y in c n\ ( &n)( ( &n\- The 
constants c n are positive and sum to 1. The density is the convex sum 
Xn c n n n( r ) or densities from each participating determinant. A similar 
expression will arise for the kinetic energy. This approach is designed to solve 
the problems of non-interacting v-representability emanating from convex 
sums of degenerate wave functions (see section XXX). 

iv. Janak's Theorem 

A very general theorem was noted by Janak [13], based on earlier work of 
Slater and Wood[14] concerning the meaning of orbital energies. Let us return 
to the functional of Eq. (7.3.11) and assume now that the occupation numbers 
n q are given and they are all non-negative and not larger than 1 and that they 
sum up to the number of electrons. Thus, for a given set of occupation 
numbers we can search for the orbitals that minimize the following 
functional: 



The equations that the orbitals must obey are still derived from Eq. (7.3.12) 
leading to the same equations as in (7.3.13), the KS equation. Now, let us ask: 
what happens to the energy when we change the occupation number of one 
of the orbitals n q by an infinitesimal amount 8n q l When we do this the "total 
number of electrons" N e changes by this amount as well. This is not a physical 
change (since electrons cannot change by non-intereger amounts) but still 
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(7.3.21) 



mathematically speaking the change can be studied. Since n q are the 
constraints and e q the Lagrange multipliers in a minimization problem we can 
use the general result of Eq. (1.11.23) that the rate at which the minimized 
function changes when the constraints change is equal to the Lagrange 
multiplier: 



dEv 
dn n 



e q (7.3.22) 



This relation, giving some meaning to the orbital energies is called Janak's 
theorem. This theorem is quite general but relies on some analytical 
assumptions of the energy functional. For example, when the occupation 
number is 1 the change can only be by a negative amount and when it is zero 
- only positive. For approximate functionals, that are analytical with respect 
to any n(r) and (p q (r) this relation holds. Such is the case for the often used 
local, semilocal and even most hybrid functionals, including Hartree-Fock 
theorem. 

In the conext of Hartree-Fock theory this result is a restatement of Koopmans' 
theorem, by which - e q is the unrelaxed ionization energy from orbital q. The 
orbital relaxation is a second order effect and thus negligible when occupation 
numbers change infinitesimally. 

E. "Virial Theorem" related identities in DFT 

The following development is inspired by the virial theorem treatment. It 
continues to consider dilation relations. Taking the derivative of n y (r) with 
respect to y, we have ^n y (r) = y 3 r • Vn(yr) + 3y 2 n(yr) and since, V • (r/(r)) = 
3/(r) + r • V/(r), we find 
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|;n y (r) = y 2 [V s -(5n(5))] s=yr (7.4.1) 

For a general density functional, A[n], we have, using the chain rule for 
derivatives, a[n](r) = SA[n]/Sn(r) and Eq. (7.4.1)then: 

d 



A[n Y ] = j a[n Y ](r)Y 2 [V s ■ (sn(s))] s=y? d 3 r = j y- 1 a[n y ](y- 1 5)V s • (sn(s))d 3 s 



= — y 1 I (r • V r a[n y ](r)) _ i n(s)d 3 s 



From which a completely general virial-dilation relation holds for any 
functional: 

^feJ + Y 2 j (r- V r a[n y ](r)) n( Y r)d 3 r = (7.4.2) 
We will be especially interested in the case of y = 1. Thus the basic relation: 

( dty ) + j(r-Va[n](r))n(r)d 3 r = (7.4.3) 



We further find that 



f, \ , / divAln] - A\n v ])\ 

A[n] + J (r • Va[n](r))n(r)d 3 r = V — 1 yJ7 j (7.4.4) 

Let us apply this for the Hartree energy. From Eq. (7.1.8) we find the following 
relation, valid for any y: yE H [n] — E H [n y ] = 0, so: 

£ H [n] + j r Vv H [n](r) n(r)d 3 r = (7.4.5) 

A similar relation, namely — i^[riy] = hold also for the exchange 
energy K x [n] (see Eq. (7.1.10)), and so: 

K x [n] + j r Vv x [n](r) n(r)d 3 r = (7.4.6) 

Where i^[n](r) = SK x [n]/5n(r) is the exchange potential. 
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Now, what about the correlation energy? Let us consider the KS DFT 
functional: 

E v [n\ = T s [n] + j v(r) n(r)d 3 r + U s [n] + E c [n] (7A.7) 

Suppose n(r) is the density minimizing minimizer of E v [n] and now plug into 
the latter the scaled density. This will give a y dependent energy: 

E(y)=E v [n Y ] (7.4.8) 

And clearly the minimum is at y = 1 so: 

= (7.4.9) 

Now, let us evaluate E'(y) using Eq. (7.1.5): 

E'{y) = 2yT s [n] + Kr)n r (r)d 3 r) + U s [n] + ^cfoj, ( 7 - 4 - 10 ) 

and plug in y = 1. We obtain: 

^[£ c [n y ] + j v(r)n r (r)d 3 r} = -2T s [n] - U s [n]. (7.4.11) 

or, using Eq. (7.4.3): 

j (r • V r (i> c [n](r) + v(r)))n(r)d 3 r = -2T s [n] - U s [n] (7.4.12) 

The second term in the parenthesis can be related to the interacting system. 
Indeed we have, using Eqs. (6.6.2): 

E[Vy] = Tf¥y] + U[Wy] 4" ^ V(f) Uy (r) ^ V 

= y 2 T[V] + yU[y] + j v(r) n Y (r)d} 

Taking the derivative with respect to y, remembering that ^"[^y] x = we 
find, for y = 1: 



(7.4.13) 
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= 2T[W] + [/[¥] +^-{\ v(r) n Y (r)d 3 r) (7.4.14) 

ay \J ) y =1 

(As a sidenote, you can see that substituting this relation in Eq. (7.4.1) gives 
after trivial manipulation the virial theorem 27^] + = j(r ■ 
Vv(r)^n(r)d 3 r). Continuing the above development, using Eqs. (7.4.11) and 
(7.4.14) we finally find: 

^E c [n Y ]\ Y=1 = 2T c + U c , (7.4.15) 

which can be written equivalently as: 

E c [n] - j-E c [n Y ]\ Y=1 = -T c [n], (7.4.16) 

or, using Eq. (7.4.3): 

W + |(r.«W>WA=-W.] S a (7.4.17) 

The Hartree and exchange energies give zero while the correlation energy 
gives a negative quantity equal exactly to —T c [ri\. 

This result shows that the correlation energy functional and the correlation 
potential are enough for determining the correlation kinetic energy (and from 
it, by T = T s + T c ) the kinetic energy itself. 

This latter result is related to the virial theorem of Slater, which shows that 
one can derive the kinetic energy of the electrons from the Born-Oppneheimer 
potential surface itself (Eq. (1.9.30)). Since DFT gives, in principle the Born- 
Oppenheimer potential surface, one can access the kinetic energy (and the 
potential energy) from the DFT calculation. 
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F . Galilean Invariance 

A basic property of the electron-electron interaction is that if the coordinates 
of all electrons are shifted by a constant L: 

r'n=r n -L (7.5.1) 

the e-e interaction enegy does not change 

1 V" 1 1 V 1 

This property is shows that the e-e interaction is translationally invariant. This 
property is also called Galilean invariance. The same property holds when the 
coordinates of each electron are rotated around some axis. This roatation can 
be described by a 3*3 orthogonal matrix 0, where T = 00 T = I: 

K = 0r n (7.5.3) 

The lengths of vectors are preserved under a rotation: 

\K\ 2 = rfrZ = (Or n yOr n = r T n O T Or n = r T n r n = \r n \ 2 (7.5.4) 

Thus the e-e interaction enegy does not change 

U " ~ 2 L \K-r'^\ ~ 2 L \0r n -0r m \ ~ 2 L \0(r n -r m )\ 

n*m n*m n*m 



2Lr„- 



(7.5.5) 



\r n -r m \ 



Thus the e-e energy is rotational invariant. 

These relations indeed hold for the density functional E H [n] since it is a 
reflection of the e-e functional. Indeed, define the shifted density: 

n'(r') = n(r) = n(r' + I) (7.5.6) 

Then: 
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1 ff n'OiKOo) „ „ 

E " ln] =2}jn^kr dr;dri 



1 ffnCr^ + LMr'2 + L) , 

— d*r[d*r 2 



-IS 



(7.5.7) 



K~r' 2 \ 

We now make a variable replacement: r„ -> r n — L and obtain: 
v r n 1 ^ ^ 73 ^3 



= 1 ff n( - r i 
2JJ In 



(7.5.8) 



■d 3 r 1 d 3 r 2 = E H [n] 



r 2 \ 

This condition, that £ H [n] = E H [n 1 ] is called Galilean invariance. It is easy to 
show that E H [n] functional is also rotational invariant. 

The exchange energy is also Galilean invariant, since the translation of the 
density will cause a translation of the density matrix: 

p{n'W lt r' 2 ) = p[nKr 1>r2 ) (7.5.9) 

It is easy based on this to show that E x [n] — E x [n']. 

Finally, the same will hold for the kinetic energies T[n], T s [n]. All this shows 
that we must demand this invariance of the correlation energy: 

E c [n") = E c [n'} = E c [n] (7.5.10) 

One consequence is the property of Galilean covariance of the potentials for 
each of the above energy functionals. For example, for the correlation energy 
we have the following result. Suppose we shift the density by a small 
displacement: 

n'(r) = n(r' + SL) = n(r') + SL ■ Vn(r') 

Thus: 
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f oE 

E c [ri] =E' c [n + 8L- Vn] = E c [n] + 8L ■ — — ^- Vn(r')d 3 r' (7.5.11) 

J on{r') 

Since we demand Galilean invariance E c [n'] — E c [n] we find: 

f SE 

51 'J ^#) Vn(r ' )dV = ° (7 - 5 - 12) 
Since ^zjk — v c [n] (r) and SL is arbitrary: 

/.JnKDVnOOd^o (7.5.13) 

For finite systemswhere the density drops to zero at infinity we can move the 
nabla sign to the potential: 

j -Vv c [n] (r) n(r)d 3 r = (7.5.14) 

— Vv c [n](r) is the force derived from the correlation potential. This shows that 
the total correlation force is always zero. Another consequence from Galilean 
invariance of the correlation energy is the Galilean covariance: 

Vc [n'](r') = v c [n] (r) = v c [n] (r' + I) (7.5.15) 

Similar conditions can be proved from rotational invariance. For example the 
torque: 

/rxVMn](r)„(r)d3r = (7-5.16) 

G. Holes and the adiabatic connection 
i. The exchange and correlation holes 

Let us now take a step back and return to wave function theory. We examine 
the electron-electron interaction energy 
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U = {% s \u\V gs ) (7.6.1) 
Which we write using the following operator: 



e V 1 e z ff n 2 (r,r') ^, , 



In the first term, we sum over the pairs of r-vectors of each of the N 
coordiantes (indices i and j). In the second term we use the definition: 

n 2 (r,r') = ^ 8(r - t{)S(r - fj), g 6 3 ) 

i*j 

which is the pair density operator. Note the relation between the one and two 
densties: 

n 2 (r,r') = / — ri)8(r' — fj) — y S(r — f £ )<5(r' — r{) 

if i (7.6.4) 

= n(r)n(r') — S(r — r')n(r) = fi(r)[n(r') — S(r — r')]. 
With n(r) = £f =0 — r d- With this definition, we have: 

j n(r)d 3 r = N 
j n 2 (r,r')d 3 r' = (JV - l)n(r) (7.6.5) 

j n 2 (r,r')d 3 r'd 3 r = N(N - l)n(r) 

The expectation value of n(r, r') is the "two electron density function": 

r(r,r') = {% s \n 2 (r,r')\W gs ) (7.6.6) 
With the two-electron density function, the interaction energy is: 

t/ = y//j^d 3 rdV (7.6.7) 
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This pair density function has the symmetry, positivity and normalization 
properties given by: 

r(r,r') = r(r',r) 
r(r, r') > 



/ 
/ 



r(r,r')dV = (iV - l)n(r) 

(7.6.8) 

r(r,r')d 3 r = (iV - l)n(r') 
r(r,r')d 3 rd 3 r' = (N - 1)JV 



r(r r') 

The normalization allows interpretation of -^j ^ as the probability density to 

find an electron at r and another electron at r' . One property that is intuitively 
expected of T(r, r') is as the limit |r — r'\ -> oo is approached electrons will 
gradually uncorrelate and T(r,r') collapse to the density product n(r 1 )n(r 2 ). 
Indeed, this is insight bears out in most cases: 

T(r, r') -> n(r)n(r') (|r — r' | -> cx>, in most cases), (7.6.9) 

but not always. For example, the ground state wavefunction of the Carbon 
atom W c (r 1 ,r 2 , —) (we are neglecting to write spin indices for sake of 
notational simplicity) in the large r x limit: for minimal energy reasons the 
remaining electrons will lower their energy to a maximal extent thus, the 

wave function should obey: x V c (r 1 ,r 2 , —) ^ l——^¥ c + ri (r 2 ,...) where 

l F c + ri (r 2 , ...) is the ground state wave function of the cation C + and I— — is 

the root of the propability to find an electron far at r 1 . Note however that this 
cation has a 3-fold degeneracy in its groundstate energy and thus for any 
finite r lt no matter how large, the l F c + ri (r 2 , ...) wave function is that 
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degenerate wave function distorts in a certain fashion in correlation with the 
direction r 1 /r 1 (for more details see Phys. Rev. A 49, 809 (1994) or J. Chem. 
Phys. 105, 2798 (1996)). 

We can also look at the conditional probability density to find an electron at r 
given that there is one at r' (this latter probability is n(r')/N e ), given by 

p(r ^ = w. - bio <7 ' 6 - 10) 

Obviously if one integrates on r one gets unity. We can thus view: 

n cond (r\r') = (JV e - l)P(r\r') = (7.6.11) 

n{j ) 

as a "conditional" density the density at r of N e — 1 electrons: all "other" 
electrons except that one electron is known to be at r'. Indeed, upon 
integration over r, we get, irrespective of r' : 

j n cond (r\r')d 3 r = (JV e - 1) (7.6.12) 

Furthermore, we have for |r — r' \ -> oo: 

n cond (r\r') -> n(r) (|r — r'| -> co, in most cases), (7.6.13) 

This shows that the density far from the localized electron is unperturbed. 
Now, let us subtract from this conditional density the total N e -electron density 
n(r) and obtain the Fermi-Coulomb hole function: 

n FC (r\r') = n cond {r\r') - n(r) (7.6.14) 

Since we localized an electron at r' the rest of the electrons will "rearrange" so 
as to be repelled from the stationary source. This will give us a "missing 
density" or "hole density", i.e. the charge density at r expelled by an electron 
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at r'. We expect the total charge of the hole is —1. Indeed plugging Eq. (7.6.14) 
into (7.6.12) we find: 

j n FC (r\r')d 3 r = -1. (7.6.15) 

Furthermore far from the hole center we have from Eqs. (7.6.13) and (7.6.14) : 

FC \ \ - -» (|r — r'l -> oo, in most cases), (7.6.16) 
n(r) 

This shows that the FC hole decays faster than the density far from the 
system: for localized systems it is, in most cases, a highly localized overall 
singly charged distribution. 

Now, because r(r.r') = n(r')n cond (r\r') = n(r')(n(r) + n FC (r\r')), the 
Coulomb interaction energy can be written as: 

U = e -\\ <r%n(r) + n FC (r\r')) 
2 JJ \r-r'\ 

And in terms of the FC hole: 

H 2 JJ \r-r'\ 
Thus the part of the interaction energy beyond the Hartree energy is the sum 

2 / I |\ 

of all interaction energies e FC (r') = — f n Fcv r w ^3 r between an electron at r' 

t> tL v / 2 \r-r'\ 

and its Fermi-Coulomb hole n FC (r\r'). We will shortly see that the correlation 
energy adimits a similar analysis only with slightly modified quantities. 

By using non-interacting electrons we can also pull out of the integral the 
exchange energy and write: 

U=U s ^\\ n V nci V^r' (7.6.19) 
5 2 JJ \r-r'\ 



Electron Density Functional Theory 
© Roi Baer 



Page 174 



This leaves a Coulomb hole which is overall neutral. It too is localized. We 
discuss this in the section after next. 

ii. The Fermi-Coulomb hole for harmonic electrons 

Let us calculate these quantities for our 2-harmonic electrons in their ground 
state triplet (so we have both exchange and correlation). The pair density and 
density for the wave function in (1.6.18) is: 



The density is plotted for several values of the correlation constant 0, 




(7.6.20) 



The density can be obtained by integrating: 




(7.6.21) 
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Figure VII-1: The 1-particle density, for a system of two harmonic fermions placed in a 
harmonic well in their triplet ground state for various interaction strengths. When 6 = - 
there is no interaction and the dip in x = is due to the "Pauli repulsion". As interaction 
grows the dip becomes deeper and broader. 

We plot the conditional probability p(x|x x ) for this system in Figure VII-2. 
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8= - fc - 




Figure VII-2: Contour plots for the conditional probability distribution p(x\x{) for a 
system of two Fermions in their triplet ground state for various interaction strengths. 
When = | there is no interaction and the only correlation is due to the Pauli principle. 
As interaction grows the probability distribution rotates by 45°. 

The XC hole is plotted next in Figure VII-3. 
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9= - e= - 




Figure VII-3: Contour plots for the FC hole n X c{x,x{) for a system of two Fermions in their 
triplet ground state for various interaction strengths. When 6 = - there is no interaction 
and the only correlation is due to the Pauli principle. 

iii. The Fermi hole in the non-interacting system 

Let us consider now the FC hole in the non-interacting system. Since there is 
no correlation in absence of interactions, we attribute the hole only to the 
exchange (Fermi) effects. A non-interacting system having the density n(r) 
that has a closed shell Kohn-Sham determinant, composed of orthonormal 
orbitals a (r), where 1 indicates spin up and -1 spin down. 
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r s (r,r') = {Hn] gs \fi(r,r')Mn] gs ) 

= ^[0a(O 2 fc (r') 2 ~ 4> a (r)(t> b (r)(t> b (r')(t> a (r')] 

a=tb 

= ^[0a(O 2 ft (r') 2 - <p a (.r)<p b (r)<p b (r')<p a (r')] 



a,b 

= n(r)n(r') — P(r,r') 2 

Where the sum is over the orbitals in 5S (the occupied KS orbitals) and we 
defined the density non-interacting matrix P(r, r') = £ a 0a( r )0a( r ')- 

Exercise: Prove that 

j P(r,r') 2 d 3 r = n(r') (7.6.23) 
Exercise: As a check, integrate V s (r, r') over r and find: 

| r s (r,r')d 3 r = n(r')[N e - 1] (7.6.24) 
The Fermi-conditional density is: 

P(r r'Y 

n(r\r') = n(r) - (7.6.25) 
n(r ) 

And the Fermi-hole is: 

P(r r'Y 

n F (r\r') = ^f- (7.6.26) 

n(r') 

It can be shown[15] that in most cases the density matrix P(r, r') decays 
exponentially as |r — r'\ -> cx>, although this could be much slower than n(r). 
Thus we may say: 

n F (r\r') -> (|r — r'| -> cxi, in most cases), (7.6.27) 

This is weaker than Eq. (7.6.16) for the total FC hole. This shows that n c 
ii r; — n F decays to zero in a similar but opposite way than n F (r\r'): 
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n F (r\r') -> — n c (r \r') (|r — r' | -> oo, in most cases), (7.6.28) 
Based on Eq. (7.6.23) the Fermi hole carries all the charge of the FC hole: 

j n F (r\r')d 3 r = -1 (7.6.29) 

This allows one to say that it is the Fermi or "exchange"-hole in the non- 
interacting system that "carries the charge" of the exchange correlation-hole 
in the interacting system. Once the interacting system has been mapped onto 
the non-interacting system the Fermi-hole is easily calculated. This can be 
used to define the Coulomb hole by: 

n c (r\r') = n FC (r\r') - n F (r\r') (7.6.30) 

It has no total charge: 

j n c (r\r')d 3 r = (7.6.31) 

The interaction energy can be written now as: 

1 ff n(r')n r (r\r') „ „ 
U = U s + - , , d 3 rd 3 r' 7.6.32 
2JJ \r-r'\ 



Exercise: Compute the Fermi-hole function of the homogeneous electron gas 
Solution: We already determined the density matrix (see Eq. (5.2.24)) 

P(s) = 3n 7l(5fcF) (7.6.33) 

SKp 
sin x x cos x 

where s = |r — r'\ and ji(x) = — . The x-hole is of course independent 

of r: 

HEG f , \ ^( S ) n A( 5 ^f) 2 /7/:a/|\ 

n F (r; r + s) = = -9n , , (7.6.34) 

n (sk F ) 2 

Plotting shows the form of the HEG x-hole function: 
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Px < s )/ n 




Here kp = 3n 2 n (unpolarized gas). Since An f" ^ 1 ^ dx = n , we find HEG 
exchange energy is, : 



Up _ 1 1 [f 
~N~~N2JJ 



U F 1 1 rfn(r)n^ EG (r;r') , , 1 f°°njF G (s) 

-: : d s rd s r =-4n 

\r-r'\ 2 J s 

'AW 2 



9n 



-47T 



?k 2 



I 



3 /3 X 1/3 



(7.6.35) 



This is indeed the LDA exchange energy per particle, 
iv. The Adiabatic Connection 

Having written down the relation between the interaction energy U and the 
XC hole, we still have no such relation for the correlation energy E c . We now 
derive such a relation. Remember the correlation energy is defined as: 



E c [n] = T[n] - T s [n] + (U[n] - U s [n]) = T c [n] + U c [n] 



(7.6.36) 



Given a ground-state electron density n(r), consider a family of iV e -"electron" 
systems, with parameter < X < 1, where: 



H x = f + j v x (r)n(r)d 3 r + AU 



(7.6.37) 



The ground-state is denoted W x . The potential v x (r) is chosen in such a 
manner that the density at the ground-state wave function is n(r), i.e.: 



m|n(r)|¥ A )=n(r) 



(7.6.38) 
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This is a generalization of the idea by Kohn and Sham, that the interacting 
electron system is mapped onto a non-interacting electron system with the 
same density. Except that now we map our system to a system of electrons 
with interaction XU. When X = we have the non-interacting system and 
Vx =Q (r) is the Kohn- Sham potential v s (r) and we have: 



E [n]=T s [ 



n] + j v (r)n 



(r)d 3 r 



(7.6.39) 



where T s [n] is the kinetic energy of non-interacting electrons. When 1 = 1 we 
have the fully interacting system and v x=1 (r) = v(r) is the actual external 
potential on the electron system and the energy is: 



(7.6.40) 



(7.6.41) 



£jn] = E v [n\ = T[n] + j v 1 (r)n(r)d 3 r + U[n] 
We can also define the obvious quantities: 

EaM = T x [n] + j v A (r)n(r)d 3 r + U A [n] 
From Eq. (7.6.36): 

E 1 -E = JbiCr) - v (r)]n(r)d 3 r + U s [n] + E c [n] 

E A -E = j[v A (r) - v Q (r)]n(r)d 3 r + XU s [n] + E C:A [n] 
With: 

EcM = (V A |f + XU\W X ) - (%|f + Af7|¥ ) 

Now, the ground-state energy of the intermediately interacting electrons 
obeys, by Hellmann-Feynman's theorem: 



(7.6.42) 



d ^ 
dX Hx 



From the second equation then: 



r) 



n(r)d 3 r + (V^U^) (7.6.43) 
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^E C:A [n] = (V A |0|V A ) - U s [n] = U C:A [n] (7.6.44) 

This expression is the differential form of the adiabatic connection. If we 
integrate it with respect to X from to 1, we find: 

E c [n] = [ Va'I^K'MA' - U s [n]. (7.6.45) 
Jo 

This formula is called the "adiabatic connection" formula for the XC energy 
[16]. We may write: E CjX = T cA + XU CiX . Then ^E cX [n] = ^T CiX + U CjX + 

a ^ u c,a and so 

dT r i\ri\ dU r i\n] 

— £^_L_L + x — £^L_L = o. (7.6.46) 
dX dX 

We can rewrite U CiX in terms of the correlation hole. Indeed, if ric(r\r'), is the 
correlation hole for the X system then using (7.6.32): 



d 3 rd 3 r' (7.6.47) 



n(r')ric(r\r') 
\r — r'| 

From which: 

1 (( n{r')n c {r\r') 
= ^ JJ | r _ r '| d 6 r'd 6 r (7.6.48) 

And we see that the correlation energy can be obtained from the the X- 
averaged Coulomb hole, called the correlation hole (since it is associated with 
the correlation energy): 

n c (r\r') = [ n£(r\r')dX (7.6.49) 

Note that because / n^(r, r')d 3 r = 0, we have also: 
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n c (r\r')d 3 r = (7.6.50) 

It is interesting that the correlation energy, like to Coulomb energy, can be 
represented as a Coulomb interaction of the density and a hole as in Eq. 
(7.6.48). Note however that the relevant hole as a coupling-constant (A) 
averaged correlation hole and not the Coulomb hole itself. 

Let us discuss one of the important consequences of Eq. (7.6.50) i.e. that the 
total charge of the correlation hole is zero for localized charge systems. If we 
rewrite the correlation energy as: 

1 ff n(r') n c (r\r') 



EM=- 2 jj 



\r — r'\ 

We see that the correlation energy can be written as 



d 3 r'd 3 r (7.6.51) 



where: 



E c [n\ = j e c (r')n(r')d 3 r' (7.6.52) 



f n ^ r 'Ks r (7.6.53) 



e< ( n c [r\V) 
ec(X) ~ 2 J \r-r'\ d 



(Note that this is just a suggestion since adding to e c (r') any function Ae c (r') 
for which j Ae c (r')n(r')d 3 r' = will give the same correlation energy). 
Because for a fix r' n c (r\r') is an oveall neutral charge density in r space its 

"Coulombic potential" r is expected to decay relatively fast for r'(faster 
than 1/r'). 

H. Derivative Discontinuity in the exchange 
correlation potential functional 
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VIII. Approximate correlation energy 
functionals 

While the correlation energy in atoms and molecules is only a small fraction 
of the total electronic energy it is found that it is in fact a very large 
percentage when one computes energy differences, such as energy of 
atomization (i.e. the difference between the energy of the atomic constituents 
to the energy of the molecule) or relative energies of different conformations, 
as those determining the shape of the Born-Oppenheimer potential surface. In 
essence, the exchange correlation energy is the chemical bonding energy. It is 
therefore crucial to model this energy accurately. We describe below some of 
the basic approximations for density functional theory. 

A. The local density approximation (LDA) 

The mapping of the interacting electron system onto the non-interacting 
system, encapsulated in Eq. XX, is of formal interest only, unless we devise a 
way to approximate the correlation potential. One way is to consider the 
correlation energy per electron e c (n) in the homogeneous electron gas of 
density n. This energy can be computed with relatively high precision using 
Monte Carlo methods. Under this approximation we can write the correlation 
energy as j e c (n(r))n(r)d 3 r. However, this does not yield in practice good 
enough results and is thus seldom used. A more successful ways was devised 
by Kohn and Sham. They considered both the exchange and correlation per 
an electron in the homogeneous gas, e xc (ri). In this case the correlation energy 
comes out: 




(8.1.1) 
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This approximation is called the local density approximation (LDA)[17]. The 
functional E X c A [ri\. It leads to the following LDA approximation for the 
energy functional: 

E{; DA [n] = T s [n] + j v(r)n(r)d 3 r + E H [n] + j e xc (n(r))n(r)d 3 r (8.1.2) 

The minimization of this functional, by the Kohn-Sham approach leads to the 
LDA approximation of DFT This approach is highly successful and is 
considered in DFT as the basis for most of the developments of other 
functionals. 

Note however that in LDA the correlation energy is extremely awkward 
looking because of the term - E x [n] . The presence of this term has detrimental 
effects which harm some of the predictions of DFT. 

i. The exchange energy per electron in the HEG 

In section XXX we discussed in some detail the Hartree-Fock theory of the 
homogeneous electron gas. We defined a Jellium as a smeared positive 
background of Volume V at density n = — together with N e electrons. We 
showed that the Jellium self energy the Jellium-electron attraction energy and 
the electron Hartree energy all cancel exactly in the HEG. Thus, the energy 
per particle is given by: 

e(n) = t s (n) + e xc (n) = t s (n) + e x (n) + e c (n) (8.1.3) 

We already calculate exchange energy, using a Hartree-Fock treatment of the 

HEG and we saw that e x = —C x n 1/3 (see Eq. (5.2.32)) with C x — -(-) E h a . 

4 \tc/ 

In terms of the Wigner-Sietz radius, which is a dimensionless quantity given 
by: 
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(8.1.4) 



is the radius of a Jellium sphere containing the charge of an electron. Thus: 



If we have a way of computing e(ri) and t(n), we can then find £xc(. n ) from 
Eq. (8.1.3). We can then also compute e c {n) from Eqs. (8.1.5). 

ii. Correlation energy of the HEG: the high density limit 

The calculation of e c can presently be done analytically in two limits. One is 
the high density limit r s -> where the kinetic energy dominates and the 
Coulomb interaction can be treated as a perturbation. In this limit the kinetic 
energy is that of non-interacting particles. Thus in the perturbative approach 

one can take H = £n=i — ~^n- The unperturbed wave function is the Slater 

Z77l e 

wave functional wave function ¥ composed of the plane wave orbitals -^=- 
with k taking the N e lowest momentum vectors k < k F . The energy of the 

ft 2 k 2 

unperturbed state is E Q = Y,k<k F The first order contribution of the the e-e 

Coulomb repulsions is (^olf/l^o). This quantity is equal to the direct and 
exchange contribution (see Eq. (4.1.12)). The direct part is nullified by the 
other electrostatic interactions, so the 1 st order contribution is essentially the 
exchange energy of the HEG which we already included e (see Eq. (8.1.3). 
Thus to continue and determine the interaction energy beyond exchange, i.e. 
the correlation energy, we must move to at least second order perturbation 
theory. When one does this, one finds that the usual second order 
perturbation theory yields infinite terms. These are associated with low wave 
length excitations where a pair of electrons having the momentum statse I k t ) 




(8.1.5) 
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and I k 2 ), are excited by the Coulomb interaction V q = 4n/q 2 to states \k 1 + q) 

1 1 

and \k 2 — q) (such that -(Z^ + q) 2 > e F and -(fc 2 — Q) 2 > ^ F ). One can show 
that for small q this process gives a term proportional to In q in the expression 
for the 2 nd order perturbation energy. This term is singular at low q. A method 
of performing perturbation theory which is non-singular and goes beyond 
second order, was devised in 1957[18]. This theory is essentially exact at the 
high density limit and leads to the following relation: 

e c = A In r s + C + 0(r s ) (8.1.6) 

Where r, = — is the Wigner-Seitz radius, namely a r s is the radius of a 

sphere in the Jelium which scoops an amount of charge equal to le, where e is 
the elecmentary quantum of charge (the electron charge). A = 0.03 llE h and 
C = -0.048 + O.OOlEft. Later work refined these constants: A = 0.031091£" h 
and C = -0.046644^. 

Exercise VIII-1 

Using Eq. (1.9.43) and the Hellman-Feynman theorem prove that 

t c = t-t s = 3v™ G - 4e™ G = 3v£ EG - 4e™ G (8.1.7) 

Hints: 

(a) Show that e 2 2 = u ee — u H where u ee and u H are, respectively, the total 
electronic repulsion energy per electron and the Hartree energy per electron. 

(b) From the the fact that e x oc n 1/3 show that the last equality is correct. 
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iii. Correlation energy of the HEG: the low density limit and the Wigner 
crystal 

The second limit is that of low density, where the electrons form a crystal. We 
give the development of this limit, originally proposed by Wigner[19] when 
he devised a theory for the electron density in metallic sodium. 

Wigner assumed that at low energy the homogeneous electron gas forms a 
crystal. Now that may sound strange: how can the density be uniform and at 
the same time the electrons form a crystal? Thanks to Quantum Mechanics 
this is actually not a contradiction, as the following example shows. 

Exercise: Show that for 2 particle in a 3D box of volume V with periodic 
boundary conditions, if the Hamiltonian is: 

H = -^l-^ 2 2+u( ri -r 2 ) (8.1.8) 
2m 2m 

Then 

2 

1) The eigenstates have a homogeneous 1-particle density n = - 

2) The pair correlation function has structure. [g{x 1 ,x 2 ) = j-_gl£3_ = j 

Pclassi x l\ x 2) n { x l) n { x 2) 

In the low density regime the electron kinetic energy (per electron) can be 
neglected since, as seen in XXX it is proportional to r~ 2 while the repulsion 
energy between the electrons per electron is proportional to r s _1 . At lkow 
density the Pauli exclusion priniciple is non-operative, since electrons do not 
overlap. Thus the quantum nature of the electron is gone at this limit and we 
can think of the electron as a classical particle that localizes. This is because 
non-localization of particles in quantum mechanics arises only from their 
need to reduce kinetic energy. The electrons will arrange themselves in the 
lowest energy state by forming a close packed crystal. Each electron is then as 
far as possible from each other electron, while still filling 3D space with 
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average density n . Let us calculate the energy of such a crystal. Consider one 
of the electrons in the crystal. We imagine it together with a cell containing 1 
unit of positive charge. This cell shape depends on the crystal symmetry. 
Following Wigner, we neglect the crystal structure and assume each electron 
is surrounded by a sphere of positive charge completely neutralizing it. Our 
approximation then neglects the volume of the space between the spheres. 
The radius of the sphere is r s and it is filled with smeared positive charge and 
with one negative charged point-electron at its center. The spheres do not 
interact since they are neutral and have no electric moments. 

The total energy per electron is the energy e S p( r s a o) to assemble the Jellium 
sphere and the energy e e i(r s a ) needed to bring the electron from infinity into 
the center of the sphere. 

Let £ sp (R) be the energy to assemble a sphere of charge density n and radius 
R. Suppose we now enlarge it by adding a shell of radius dR. The electric 
potential at distance r > R outside the sphere is Q/r where Q = ^-enR 3 is the 

charge in the sphere. The charge in the shell is dq = n4nR 2 dR and bringing it 
from infinity where the potential is zero to its place on top of the existing 

sphere involves the energy de sp = dq^ — ^—^—R A e 2 n 2 dR . Thus, by integration 
from to R, we find: e sp (R) = (~~) ~ R 5 - And so at R = r s a : 

e sp (r s ) = ~E h (8.1.9) 

Next, we want to calculate the energy to bring an electron from infinity to the 
center of the sphere. This will be done in two stages, first bringing the electron 
from infinity to the rim of the sphere, a distance a r s from its center and then 
from the rim to the center. Accordingly write e e = e rim + e center- The first part 
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is easy since we already know the potential, and it is negative since energy is 
released by this process, so: 

e rim = -^E h (8.1.10) 

Inside the sphere, at a distance R from the center there exists an electric field 
due to the Jellium, which according to Gauss' is E = r3 = —nR- This force 
is a Harmonic force, with force constant k H = — n = - 1 — . The work to move 

3 (a r s )3 

1 

an electron in this field to the center is: ^center = 

/ E(R) ■ dR = - —. The 

r s 2r s 

energy for the second stage is therefore e e i(r s ) = — -— and the total energy per 

2 TV 



electron in the crystal is: 



This then is the exchange-correlation energy for low density. We neglected the 
volume between the spheres. The exchange energy we already know from 

45817 

(8.1.5), is £x(Ts) = — : E h . Thus, Wigner's approximation for the 

correlation energy in the low density limit is: 

0.44183 

e C (r s ) = e(r s ) - e x (r s ) = E h (low density, r s oo) (8.1.12) 

r s 

Wigner also considered the correction due to the finite kinetic energy when r s 

is finite. Since we saw that the electron inside the spherical Jellium drop is a 

Harmonic potential, one can reduce the correlation energy by the 3D 

3 n — 3 



Harmonic zero-point potential, -y k H = — ^ • The correlation energy is then: 

2 2 T - 



0.44183 3 

e c (Xs) = En + (low density, r s oo) (8.1.13) 



2r; 
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iv. Monte-Carlo determination of the correlation energy for the HEG 

Between these high and low limits there is no analytical theory in general and 
a numerical computation can be made based on quantum Monte Carlo 
methods. The results of the calculation are then fitted to an analytical form 
which respects the limits 

v. The polarized HEG; local spin-density approximation (LSDA) 

Up to now we have assumed that the electron gas is unpolarized, i.e. the total 

z component of spin Sz per electron is zero. However, Sz is a good quantum 

1 1 

number and it can vary continuously from — - to -. The extreme case is the 

fully polarized case. In general one may define the density of spin-up electron 
n^(r) and that of spin down n^(r). Then: 

n(r) = n T (r) +n l (r) 

n,{r)- ni {r) ( 8 - L14 ) 

c(r) — 

n T (r) +7^0) 

For a fully polarized gas <" = 1 and the difference is first of all in the Fermi 
energy. For the HEG, since every momentum state can populate only one 
electron, we find by a similar analysis as in the unpolarized case: 

k FU = k F {i ± o 1/3 

_ T ^ _ 3 

The total kinetic energy us the sum of up and down contributions: t = 
n T t T +nj,t^ w j 1 - c ] 1 - g evaluated to be t = - EEEllfgiB us i n g the expressions for n n 

we obtain 



Electron Density Functional Theory 
© Roi Baer 



Page 192 



As for exchange energy, since exchange interaction occurs only between like 
spins, we XXXXX 



vi. Successes and failures of LSDA 

vii. Plausible reasons for the success of LSDA 

One of the uses of Eq. (7.6.48) expression is to explain the success of a simple 
theory such as LDA[20]. TO see this, let us expand the XC hole in terms of 
moments around r: 



oo 




(8.1.17) 



1=0 m=-l 



where: R = Ru = (i? sin 9 cos , R sin 9 sin , R cos 9) and 



ni m (r;R)=\ I n$ c (r,r + Ru) Y lm (9 , cp) sin 9 d9 d<p 



(8.1.18) 



Then consider the XC energy, it can be written as: 




(8.1.19) 




Were: 




o 



(8.1.20) 



And: 
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V4lr n 00 (r,R)R 2 dR = -1 



(8.1.21) 



o 



Thus, only the 00 moment of the X dependent XC holes enters the expression. 
Therefore, in a sense the angular shape of the average XC hole gets averaged 
over and only the radial dependence affects the XC energy. This is used to 
explain some of the success of LDA. 

We see that the XC anisotropy of the XC hole around r is averaged over. Only 
this average enters the XC energy formula. In LDA we use the homogeneous 
electron gas to compute the HEG. Of course this leads to an isotropic XC hole. 
Yet, since only the spherical average of the hole enters into the XC energy, this 
drastic approximation gives a reasonably good XC energy. 

Exercise: Calculate the spherically averaged X-hole for a 1-electron system (H 
atom for example) 



Solution: The orbital is ip(r), the density is n(r) = ip(r) 2 and the DM is 
p(r,r') = xp(r)xp(r') thus: 



The hole is independent of the reference point r. 

One of the important results shown below is that only the spherically 
averaged hole enters the XC energy. Thus we only need the spherically 
averaged hole: 



n x (r; r') = — n(r') 



(8.1.22) 




(8.1.23) 



Which becomes: 




(8.1.24) 
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For the H-atom n(r) = A 2 e ar , defining B = a^r 2 + s 2 y = - J s ' therefore 

v J ° r ' (r 2 +s 2 ) 

B-yJl ± y = a\r ± s\ and yB 2 = a 2 2rs : 



n s x Alel (r;s) =^A 2 j e-" |r+s| dn s = ^A 2 j e - a ^ r2+s2+2rscos9 d cos 9 



A 2 rv 
= _ e -Pji+y 



dy 



(8.1.25) 



A 2 



[(1 + a\r-s\)e~ al 



a 2 2rs 

-(l + a|r + 5|)e- a|r+s| ] 

The form of this spherically averaged hole function as function of s and r is 
shown here: 




The hole has a cusp at the r-s origin. 
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B. Semilocal functionals and the generalized 
gra di en t approxima tion 

IX. Generalized Kohn-Sham approaches 

A. The generalized Kohn-Sham framework 

The Kohn-Sham approach is based on pure density functionals for the 
exchange correlation energy. This approach is based on the Hohenberg-Kohn 
universal density functional F HK [n] is written as: 

F HK [n]=T[n] + U[n] (9.1.1) 

We noted one problem with F HK is that it is not defined for all densities. This 
problem can be cured by the Levy-Lieb procedure, where the F HK is replaced 
by a functional based on a constrained minimum principle: 

F LL [n] = min(¥|f + 0\w) (9.1.2) 

Similarly for the non-interacting system, one finds: 

T s [n] = min(o|f |o) = min5[cl)] (9.1.3) 

cj)->n 4>->n v ' 

Where here the search is over all determinants. Here we defined here a new 
determinantal functional S[<&] with the intent of generalizing the KS 
approach, as we do now. Using the minimizing determinant O s , the Kohn- 
Sham correlation energy: 

E c [n]=F LL -T s -U s (9.1.4) 

Is written using T s = (O s |f |O s ) and U s = (O s |f?|O s ). 

Since a Slater wave function is given in terms of orthonormal one-electron 
orbitals <p n (x), it is best to view S as a functional of the orbitals and we write 
S [{</>}]. We can now generalize and demand that S[{0}] not be just the kinetic 
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energy but a more general functional of the Slater wave function. We denote 
the functional derivative of S as follows: 

SS [<D] 

= O s n (r) (9.1.5) 



8<Pn(r) 

Where S is some convenient operator. 

We now define the energy functional for a given potential v s (r) and 
determinant O: 

E s [v Sl O] = 5[0] + J u s (r)n[0](r)d 3 r (9.1.6) 

Where: 

n[0](r) = (0|n(r)|0) (9.1.7) 

We then minimize with respect to the orbitals of 

O, keeping them normalized. This constrained minimization results in the 
orbital equations: 

(ds + v s (r))cp n = e n( p n . (9.1.8) 

Furthermore, let us define the density functional F s resulting from 
minimizing S with constrained of a given density: 

F s [n] = minS[<D]. (9.1.9) 

From this we further define the density functional i?s[n]: 

F LL [n] = F s [n] + R s [n] (9.1.10) 

In actual implementations the functional R s must of course be approximated. 
Now we put all ingredients together. The HK DFT energy is developed as a 
train of equalities as follows: 
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E v = m\n^F LL [n\ + j v(r)n(r)d 3 r^ 

= mm^F s [n] + R s [n] + I v(r)n(r)d 3 r^ 

n ^ N (9.1.11) 
= min]min5[0] + fts[ n ] + I v(r)n(r)d 3 r\ 

n->N (.<I>-m J ) 

= min|5[d>] +i? s [n[d>]] + j v(r)n[0](r)d 3 rj 

In the last step we changed the order of the minimum procedure, assuming 
this is OK. Note that the final procedure is minimization with respect to 
orbitals instead of to density. The minimization, under normality of the 
orbitals (p n (x) from which the Slater wave function is built is gives the 
generalized Kohn-Sham equation: 

(6 S + v R {r) + v(rj) 4> n = e n cp n (9.1.12) 

Where: 

If we write: 

v s (r) = v R (r) + v(r) (9.1.14) 

we find that the orbitals also obey equations (9.1.8). Thus, the density n(r) 
obtained from the orbitals of the GKS procedure is the minimizing density of 
the following density functional 

E s [v s ] = mm[F s [n] + f v s (r)n(r)d 3 r] (9.1.15) 

The strength of the GKS scheme is that it incorporates explicit orbital 
functionals. This greatly enlarges the scope of approximations for DFT 
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B. Kohn-Sham from generalized KS 

By choosing: = (o|f + f/|o) one immediately obtained the usual KS 

theory with R s [n] = E c [n]. 

C. The Hybrid functional of Becke 

Let us return to the adiabatic connection formula Eq. (7.6.45): 

EM = ( {W x \OfP x )dA - U s [n] (9.2.1) 

Becke [21], suggested the following trapezoidal approximation for the 
integral: 



Jo 



^WWjWW (9 . 2 . 2) 



2 



Leading to the Becke's Half & Half approximation: 

E c in] - \ Ul ~ ~ Us = ~ CUh + Eh) ~ \k (9-2.3) 
The function (JJl + E H )[n] is next approximated as a local density functional: 
(t/c 1 + £ H )[n] « | u xc (n(r))n(r)d 3 r (9.2.4) 

Where u^ c (n) = e xc (n) — t(n) + t s (n). the interaction energy per electron in a 
homogeneous electron gas is u xc {n). This function can be computed using the 
Quantum Monte-Carlo results concerning the XC energy. 

Becke showed that this half and half theory gives a big improvement over 
LDA and in some cases also over GGA functionals. This is because he has 
found a way to approximate Hartree-Fock theory and LDA. 

The Becke Half & Half orbital (BH&H) was the first example of a hybrid 
theory a theory which is a mixture of a Hartree-Fock functional of orbitals 
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and a DFT functional of density. Becke eventually wrote down a functional 
which contains 3 parameters [22]: 

E B3LYP [0H -'^N e ] 

* E#p A [n] + a (ff - Ek SDA [n]) + %(f| 88 - f£ SDA [n]) (9.2.5) 
+ a c (E p c ^[n]-Eh SDA [n]) 

We see that the B3LYP functional starts from a LSDA functional and adds to it 
3 corrections. One is a fraction of the explicit orbital exchange correction 
(K - £x SDA [n]) then a correction (£| 88 - E^ SDA [n]) to the LSDA exchange 
given by a GGA type functional called Becke88 [23] and finally a GGA 
correction for the correlation energy given by the PW91 functional [24]. The 
values of the 3 parameters a = 0.20, a x = 0.72 and a c = 0.81 were found by 
optimizing the performance of the functional for 56 atomization energies, 42 
ionization potentials (calculated by ASCF method) and 8 proton affinities. 



E . Long-Range self-repulsion and lack of 

derivative discontinuity 

F. Range separated hybrids 

G. Orbital functionals and optimized effective 

potentials 

H . Approximate correlation functionals and the 

Born-Oppenheimer force on nuclei 

The Kohn-Sham density functional theory method (KS-DFT), when applied 
for electrons of molecules makes sense only for frozen nuclei, because the 
external potential in DFT is assumed time-independent. However, the nuclei 
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in molecules are not in general motionless. How do we KS-DFT a useful 
approach for molecules? Here we use the Born-Oppneheimer approximation, 
which allows us to divide the electron-nuclear problem into 2 stages: first 
computing the electronic energy E({R}) for each nuclear configuration 
{R} and then combining it with the nuclear repulsion energy V rep ({R}) to 
obtain an effective potential for the motion of the nuclei. The nuclear-electron 
potential is given by: 



In order to compute the electronic energy we first write down the functional: 



We then minimize it under the constraint for the number of electrons 
J n(r)d 3 r = N. This gives a minimizing density n*(r)and we have: 



E[{R}] = E vim [n t ] = T s [n*] + v nuc (r, {R})nXr)d 3 r + E HXC [n t ] (9.3.3) 



Since in exact DFT £"[{#}] is the exact electronic energy we are assured that: 



Since this relation is true for the electronic Hamiltonian (as can be seen from 
Hellman-Feynamn theorem). 

However, what happens when we use an approximation for E HXC , as done in 
all applications of KSDFT? Does this relation still hold? We now show that 
indeed it does. 

The density n* is determined by: 
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(9.3.1) 




(9.3.2) 




(9.3.4) 



SE vm) [n] 
Sn(r) 



The electronic force on the nuclei is given by: 
F I = -V,£ = -ViE v(m [n t ] 

" * l^nuc 

(r,{R})n*(r)d 3 r 



-I 
-I 



n] 



Sn(r) 



Vjn*(r)d 3 r 



n* 



The second integral is zero because: 
8E vm) [n] 



Sn(r) 



Thus we find that the relation: 



7 i = j -V,v nuc (r,{R})n*(r)d 3 r 



(9.3.5) 



(9.3.6) 



V l n t (r)d 3 r = j [iV l n*(r)d 3 r = j n t (r)d 3 r = (9.3.7) 



(9.3.8) 



still holds. Even if the approximate XC functional does not yield the exact 
density the formal relation between the force and the density is still valid. 
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X. More on the DFT correlation energy 



A. Approximations to E x are approximation to E c 

In Kohn-Sham DFT the ground state energy of a system of electrons having 
the particle-density n(r)is written as a sum of density functionals[17]: 

E = T s [n] + f v(r)n(r)d 3 r + E H [n] + E xc [n] (10.1.1) 

Here, T s [n] is the kinetic energy of a system of non-interacting fermions having 
the same mass and ground-state density as the electrons, v(r) the external 

potential exerted on the electrons, E H = ^jj d 3 rd 3 r' is the so-called 

Hartree energy of the system and E xc [n] is the XC energy functional. The 
explicit form of the latter functional is formally known but practically 
inaccessible to computation. The basic challenge in DFT is to develop an 
efficacious approximation to the XC energy functional. 

The most widely used functionals of DFT fall into two broad categories. 
Semilocal functionals which are based locally on the density the orbitals of the 
non-interacting system or on their derivatives, such as the local density 
approximation (LDA), generalized gradients approximation (GGA) and Meta- 
GGAs [23, 25] The second category includes non-local functionals such as the 
widely used B3LYP[22]. Recently a new brand of hybrid functionals was 
developed, designed to remove the long range self-interaction in DFT[26-28]. 

For a given density one can use the wave function ^¥ s of the non-interacting 
system and define the "exact exchange energy" E x [n] = (^slf/l^s) — E H/ 

1 1 

where U = -Y,n*m — is the repulsion energy operator. The functional defined 
by E c [n] = E xc [n] — E x [n] is called the correlation energy functional. Kohn 
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and Sham[17] recognized that in practical calculations it is beneficial to 
approximate the exchange energy than to use its exact form. Indeed it was 
found that the errors introduced into E xc by using simple approximations to 
E c [n] are partially cancelled by the approximate treatment of the exchange 
energy. This error cancellation is the one of the major reasons for success of 
DFT. Since E x is actually known exactly and easy to calculate, it is evident that 
using an approximate exchange functional E^ pprox j n place of the exact £V is a 
way to correct the approximate correlation energy . For example, the 
correlation energy functional in the LDA is NOT the LDA of the correlation 
energy E^ EG = f e^ EG {n{r))n{r)d 3 r, but instead: 



E L r DA = E, 



HEG 



+ 



j ejj EG (n(r))n 



(r)d 3 r — E x 



(10.1.2) 



This view emphasizes the fact that LDA of local exchange is a way of 
correcting the correlation energy for inhomogeneous systems. Similar 
expressions are to be used in more advanced functionals (GGA etc). Any 
correction to e x for example, using a GGA type expression 
/ /x(n(r), |Vn(r)|) d 3 r should be viewed as an essentially correlation 
correction. 

This way of thinking leads naturally to hybrid functionals. A hybrid 
functional is a further improvement of LDA (or GGA) obtained simply by 
multiplying the correction in (10.1.2) by a factor A: 

e hyb,i = j ^( n ( r ))n(r)d 3 r + l J e x ,EG (n(r))n(r)d 3 r - E x (10.1.3) 
Or, based on BLYP, one of the GGA functionals we can write: 



,HYB,2 



j f c LYP (n(r),\Vn(r)\)d 3 r + X j //(n(r), |Vn(r)|)d 3 r - E- 



(10.1.4) 
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Becke [22] found by a series of calculations on many molecules that the 
"natural" assumption X = 1 is not optimal and X should be somewhat smaller, 
around 0.8. This lead to development of the extremely successful B3LYP 
functional. However, experience has shown that some systems require lower 
values of X while others strive for higher values. An example is reaction 
barriers where lower values of X (around 0.5-0.7) were needed. There is 
presently no general systematic way to set the value of X. 

The hybrid DFT approaches are well established electronic structure methods 
however in the past few years it has become increasingly clear that several 
types of calculations, for example, polarizability, electron transfer excitations, 
Rydberg states and charge allocation in weakly coupled systems are often not 
well described by the above DFT/TDDFT methods[29-31]. This deficiency, 
which is not cured by the Becke approach to hybrid functionals, was 
attributed to spurious self interaction[32] and missing derivative 
discontinuities[33] - two pervasive problems in density functional theory 
(DFT) that are intimately related. [34, 35] 

One way to mitigate the spurious self interaction and to retain a good 
treatment of correlation is to deploy a range-separated hybrid functional[26- 
28, 34, 36] In this approach, the exchange term in the Kohn Sham energy 
functional is split into long-range and short-range terms, e.g., via r _1 = 
r~ x erf(yr)+r~ x erfc(yr). The short-range exchange is represented by a 
local potential derived from a Semilocal functional. The long-range part is 
treated via an "explicit" or "exact" exchange term. If one assumes that an 
appropriate choice for y is system independent, its value can be optimized 
using a molecular training set for optimizing its value. Such semi-empirical 
approaches, typically with y in the range of 0.3-0.5 ao' 1 , were shown to 
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achieve impressive results for the ground state properties of some classes of 
systems. [27, 37, 38] Furthermore, it was demonstrated on the benchmark 
model of Dreuw and Head-Gordon[39] - the C2H4 C2F4 dimer at large 
molecular distances - that the range-separated hybrid corrects the principal 
deficiencies of the charge transfer excitation prediction of TDDFT.[38, 40] 

One of the drawbacks of the semilocal DFT and Becke schemes is that the 
correlation energy of Eqs. (10.1.3) - (10.1.4) is much too nonlocal due to the 
presence of the nonlocal exchange E x . This is unbalanced in view of the 
presence of E H in the total energy. This is the root cause of long-range "self 
repulsion" which exists in all these functionals. A way to improve the 
situation is to include in the exchange functional a non-local long-range 
exchange, which leads to the following form for the correlation[26, 27, 41]: 

E HYB,d = j ^P (n(r) | Vn(r) | )d 3 r 

+ f (//(n(r), |Vn(r)|) - / x , y (n(r), |Vn(r)|)) d 3 r ( 10 - L5 ) 

- - E X y) 

The functional E Xy is the exchange energy of a system of particles at density n 
mutually interacting through the repulsive potential u y (r) = — - — where y is 
a parameter: 

E XiY = ~\\\ \p(r,r')\ 2 u Y (r)d 3 rd 3 r'. (10.1.6) 

f x — f Xv is a semilocal density approximation to the complementary 
interaction y y (r) = - — u y (r) = — - — . The idea behind this appealing 
approach is to eliminate in the correlation functional the long range 
dependence on the density, y is treated as a universal parameter determined 
by fitting the DFT predictions based on Eq. (10.1.5) to a large database of 
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experimental data. This forms the basis of several recent semi-empirical 
approaches all of which use a universal long-range parameter y determined 
by semiempirical fitting to known experimental and sometimes high quality 
ab-initio results [26-28, 34, 36-38]. These approaches have been quite 
successful for a broad variety of molecules, sometimes beyond the types used 
in the fitting procedures. 

The approach of a universal y exchange functional with complementary 
semilocal functionals cannot however claim universality. In some systems the 
exact value of y is critical and no semilocal DFT method can correct for it. We 
show examples below. An alternative view was developed in ref. [28], which 
provides a theory for the long-range parameter y based on the adiabatic 
connection theorem, from which the following exact expression for the 
correlation energy can be deduced: 

E c = (V\? Y \V) - (V s |7 y |V s ) (10.1.7) 

Where Y Y = -Y, n ^rny Y (Xnm)- I n rer - [28] arguments were given for the 
existence of a parameter y for which Eq. (10.1.7) holds for practically all 
densities. This equation describes the correlation energy as a difference 
between potential energies of interaction in the interacting and in the non- 
interacting systems. Unlike the standard definition of correlation energy as 
E c = (v\f \v) - (^sl^l^s) + - (^sl^l^s) - Ex, the there are no 

explicit contributions of kinetic energy differences in Eq. (10.1.7). These are 
absorbed by the interaction screening parameter y. In principle all quantities 
in this equation, including y, are density-dependent. We stress that the 
difference between this exact theory and the usual exposition of RSHs is that 
the latter do not derive (or even strive to derive) an exact expression for the 
correlation energy. Eq. (10.1.7) can be applied to the HEG and the exact value 
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of y can be determined as a function of density [38, 42]. The result is shown in 
Figure X-1 and it is seen that 7 is in fact strongly density dependent, certainly 
not a universal constant. In particular, for molecular densities, where r s is 
between 0.5 and 5 the long-range parameter 7 changes by an order of 
magnitude. 

When we add the exact exchange 
energy to the exact correlation 
energy of Eq (10.1.7) we obtain for 
the exchange-correlation the 
following expression[28]: 

Exc 

= (¥|f y |¥) 

--j n{r)n(r')y Y {\r (10.1.8) 
-r'\)d 3 rd 3 r' + E Xy 

— pY 1 p 

— n xc " r n X,Y 

This expression is exact. The right- 
hand equivalence defines a new functional, the complementary XC energy 
functional E y xc . It is the focus for approximation in our proposal. Because the 
potential y Y (r) is short-range, we expect that a semilocal approximation, 
depending on density and perhaps gradients, is of sufficient accuracy for 
describing the "complementary" XC energy. 

Except for the HEG, we have no practical way of calculating the expectation 
values in Eq. (10.1.8) so an approximation must be developed for determining 
the long-range parameter 7 and for determining the complementary XC 
energy functional EZ C . The procedure for determining the variable 7 will be 
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Figure X-1: The parameter y as a function of 

the density paramter r s = a It-w) for 
the fully Polarized and unpolarized 
homogeneous electron gas. 
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called below "long-range parameter tuning". We now show some preliminary 
results which demonstrate the fact that one cannot rely on a universal long- 
range parameter y. 

B. The necessity of long-range parameter tuning 

The requirement for a system-dependent parameter is not just a theory it 
comes up in practical calculations. Here we give two examples of unpublished 
results, one in DFT and one in TDDFT. 

i. The symmetric radical cation 

In DFT, a clear example of this issue is seen in the symmetric radical cation 
systems of the type R% -* R + R + , where R is any molecule or atom. Specific 
examples are R = H, He and Ne. It is known that for these systems a semilocal 
or even nonlocal (B3LYP) DFT gives qualitatively spurious potential surfaces 
and the big issue is self-repulsion[l, 43] .Yet, even with a standard y-functional 
the results are not satisfactory. Take for example the bond dissociation energy 
(BDE). This can be estimated in two ways. The first is the atomization energy. 
The second is the depth of the molecular well relative to the asymptote of the 
potential surface. From both numbers we subtract the zero-point vibrational 
energy to get the BDE. We can see in Table 1 that the BDE of R=He and Ne 
based on atomization energy using a semi-empirical y = 0.5 1 (denoted 
BNL functional[38]) gives BDE around 1.5-2 times the experimental value! 
While in the asymptotic method the calculated BDE is too small by 15%. Such 
discrepancies cannot be waived by a better semilocal functional. They reflect 
an inherent imbalance between short and long range exchange inflected by a 
wrong long-range parameter y. Indeed, by tuning the long-range parameter in 
the following way we obtain a value y* which is specific for each system. The 
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tuning procedure is "ab initio" in that it does not use any experimental data, 
and results from the calculations themselves. It is describe in ref. [1] and is 
equivalent to the 
procedure depicted 
in a more general 
treatment in ref 
[38]. 



Table 1: (Taken from ref [1]) Data for R 2 + , R=H, He and Ne 
calculated in the cc-pVTZ basis[2, 3]. All energies in kcal/mole. 
r eq is the inter-nuclear distance at which the E(r) is minimal. 



oi = ^2fi 1 E"(r eq ) is the harmonic frequency and "BDE" is 

Eatoms ~ E( r eq) + l h <*> where E atoms is estimated either as the 
sum of atomic energies (SUM) or by the asymptotic value of 
the potential energy curve (ASS). a R is the calculated 



-E'(r)r 5 is the 

2 v J 

polarizability as estimated from the form of the asymptotic 



polarizability of the atom R and a eff = Iim,_, 



potential. Finally a d ^ UB is the polarizability of atom R 
computed with a large basis-set. 



The resulting 
functional is 
denoted BNL*. This 
procedure chooses 
y* in such a way as 
to impose the 
"energy removal 
theorem" namely 
that the energy to 
remove an electron 
from the interacting 
system is equal exactly to —e H0M0 (the HOMO energy in the non-interacting 
system). [44] The parameter y* , as we compute it, obeys the equation 
y* — J7Y* Yivn _ cyVw _ ^^ r y 



Property 


R 


BLYP 


B3LYP 


HF 


BNL 


BNL* 


Exp. [4] 


BDE by E\' 


H 


66 


65 


60.9 


60.9 


60.9 


61 


He 


82 


75 


43 


74 


59 


55 


J atoms 


Ne 


75 


60 


2 


59 


34 


32 


BDE by E\ ' 

J atoms 


H 


NA 


NA 


60.9 


50 


60.9 


61 


He 


NA 


NA 


43 


42 


59 


55 


Ne 


NA 


NA 


2 


27 


34 


32 




H 


1.1 


1.1 


1.06 


1.2 


1.06 


1.05 


He 


1.2 


1.1 


1.075 


1.2 


1.078 


1.080 


Ne 


1.9 


1.9 


1.7 


1.760 


1.72 


1.765 




H 


2.7 


2.9 


3.3 


2.9 


3.3 


3.32 


He 


1.7 


2.0 


2.5 


2.1 


2.5 


2.42 




Ne 


0.5 


0.6 


0.9 


0.726 


0.8 


0.729 




H 


NA 


NA 


1 


0.6 


1 


1 


a n/% 


He 


NA 


NA 


0.98 


NA 


0.98 


1 


Ne 


NA 


NA 


1.01 


NA 


1.02 


1 


<"'{<) 


H 


5.3 


5.6 


4.51 


5.8 


4.51 


4.50 


He 


1.6 


1.5 


1.34 


1.8 


1.41 


1.38 


Ne 


3.1 


2.9 


2.4 


3.2 


2.70 


2.66 



- 6 



H0M0 = EY\N)-EY\N - 1) where e y H0M0 is the HOMO energy and E^is 



the calculated ground-state energy of the N -electron system. This procedure 
approximately imposes the energy removal theorem and defines BNL*. The 
merit is that the HOMO and the SCF procedures are now well-balanced. 
Unlike the results of BNL, where y is not optimal, those of BNL* are very 
satisfactory for several observables when compared to experiment, as seen in 
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Table 1. The comparison with regular functionals is convincing and the 
improvement over BNL and Hartree-Fock is also large. Even the complete 
basis set limit of the atomic polarizability as calculated by BNL* compares 
accurately with experiment. 

ii. Aromatic donor -TCNE acceptor charge transfer excitation 

The second example concerns charge-transfer excitations in TDDFT. The 
results shown here are part of a manuscript recently submitted for 
publication[45]. It is known that semilocal functionals cannot reasonably 
describe such excitations[39], so a range-separated hybrid is a welcome cure. 
We tested our approach on complexes formed by an aromatic donor 
(Ar=benzene, toluene, o-xylene and naphthalene) and the tetracyanoethylene 
(TCNE) acceptor, for which optical absorption is available both in gas phase 
and in solution. [46] All calculations were performed using QCHEM 3.1, [47] 
modified to include the range-separated BNL functional, [38] using the cc- 
pVDZ basis set. [3] The internal structure of the molecules in the complex is 
known to be little-perturbed by complex formation[30, 31, 48, 49] and the 
equilibrium distance and relative orientation of the 7t-stacked donor and 
acceptor determined from the conventional B3LYP[21] hybrid functional is 
known to compare well with experiment (where available). [30, 49] Therefore, 
B3LYP-optimized geometry was used throughout. The TDB3LYP excitations 
energies were much too low and the BNL much too high (see Table 1). In 
molecular complexes, the lowest photon energy required to induce a CT 
excitation, h vct, is given for asymptotically large donor-acceptor distances by 
the Mulliken rule [50] 

hv CT = IP{D) - EA{A) - 1/R (10.1.9) 
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Where IP(D) and EA(A) are the donor ionization potential and acceptor 
electron affinity, respectively. The last term on the right hand side is the 
Coulomb energy of attraction between the electron-hole pair formed by the 
charge transfer, where R is the inter-molecular separation. For our TDDFT 
calculation to conform to the Mulliken rule, the ionization energies computed 
from Eq. (1) must correspond to the HOMO energy of the neutral donor, but 
also to the HOMO energy of acceptor anion. Thus, one needs to generalize Eq. 
(1) so as to yield, as closely as possible, two limits. We therefore look for y that 
minimizes the following / (y) function 

/G0= Yu H'omo ~ {EyiN, - 1) - E Y m)\ (10.1.10) 

i=D°,A~ 

For complexes where a range-parameter y that renders / (y) very small can be 
found (which is the case for all complexes we studied), we expect the range- 
separated hybrid to yield a quantitative description of CT excitations. An 
example of the way this tuning procedure works and its ability to replicate the 
Mulliken law for large separations R in sown in Figure X-2 




0.2 0.3 0.4 0.5 0.01 0.02 0.03 0.04 0.05 

YCao 1 ) 1/R(a0) 

Figure X-2: Tuning y in the Benzene-TCNE complex. Left: / of Eq. 
Error! Reference source not found, as a function of y. The optimal y is 0.331 for which / is 
very small. Right: The Mulliken rule, compared to TDDFT results obtained from the 
optimal y, as a function of inverse inter-molecular separation. 
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With the asymptotic behavior enforced, we expect proper balance between the 
semilocal and non-local exchange components. In Error! Reference source not 
found, we compare the calculated and experimental gas-phase results for 
various Ar-TCNE complexes. It is readily seen that the B3LYP results are 
unacceptably low and predictive power is absent. Results of generalized- 
gradient (GGA) calculations (not shown) are even lower than the B3LYP ones. 
With the range-separated BNL functional, [38] excitation energies determined 
with an "off-the-shelf" y value of 0.5 ao' 1 are much too high with respect to 
experiment. But with y* quantitative agreement is obtained to within ±0.2eV. 
For benzene and toluene, the theoretical oscillator strengths are also in good 
agreement with experiment, but are too weak for xylene and naphthalene. 
This is likely a basis-set issue, as oscillator strengths are much more sensitive 
to the basis set than the excitation energies. 
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XI . TDDFT 

A. Time-dependent Linear response theory 

The linear response of 

ih—xp k (r,t) = (fi\p(t)] + v 1 (r,tj)ii} k (r,t) 

dt V J (10.2.1) 

Where the TD Hamiltonian, in the adiabatic TDDFT approximation, is: 

H[p] = f + v(r) + v H [n(t)](r) + i# c [n(t)](r) + RY[p] (10.2.2) 
And the TD density matrix is defined in terms of the TD orbitals: 

N 

p{r,r',t) = Y j ^(T,t)xp k (r',ty (10.2.3) 

k=l 

The initial orbitals are assumed the lowest eigenstates of H[p ]: 

8\p°M(r) = eM(r) (10.2.4) 

where: 

N 

p {r,r')=Y j r k {r)r k {r'T (10-2.5) 

k=l 

We assume that the perturbation v x {r, t) is weak and we will treat it as a first 
order quantity. Before we do this however, let us prepare the equation for 
linearization by defining: 

Mr.t) = e-W* (l>Ur) +^(r,t)) (10.2.6) 

and: 
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p(r,r',t) = p°(r,r') + p\r,r',t) 
n(r, t) = n°(r) + -n}(r, t) 
With these definitions we have: 



(10.2.7) 



ft-[e-^/»(^(r)+^(r,t))] 

= (//[p°+p 1 (t)] ( 10 - 2 - 8 ) 
+ v\r, 0) [e-fctf/ft ( ^°(r) + 0£(r, t))] 

From this we find: 

ihy^lir, t) = (H[p° + p\t)] -e k + v\r, t)) (i/-°(r) + ^O, 0) ( 10 - 2 - 9 ) 

This equation is rigorously equivalent to the first equation. It is a non- 
homoegenous linear equation for the unknown ip k (t). 

We now make the linearization step. We neglect all quantities beyond first 
order. As a first step we have: 

ih-ipl(r,t) « (tf 1 + vHr.t)) Vk(r) + (fi[p°] - e k )^(r,t) (10.2.10) 

Here, the symbol "~" means "equal to first order". We use the definition: 

H 1 « tf[p° + p^t)] - H[p°] (10.2.11) 

An efficient method to obtain linear response of a perturbation is as follows. 
Select a small parameter e and solve the following time dependent equation 
(using a standard propagator, such as Runge Kutta)[51]: 



ih — xpl(r,t) = I + v 1 (r,t) \ ipk(r) 

+ {n\p ]-e k )xl>l(r,t) 

With: 



(10.2.12) 
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N N 



p\r,r', t) * £ xp° k (r)xpl(r' , t)* + £ 0j(r, t)^°(r') 

/e=l /e=l 

nHr, t) « 2fle ^ t/;°(r)^(r, t) 

/c=l 

The size of £ must be so small that the orbital amplitudes are linear with e. 
The solution involves repeated applications of the Hamiltonian. This could be 
considerably more efficient than the algebraic method described below, 
requiring the computation of all the eigenvectors of the Hamiltonian. Of 
course, for small systems with a small basis the algebraic approach may be 
extremely effective thus we discuss it now. 

B. Algebraic Approach to Linear response 

In order to derive an algebraic approach, we need expressions for all the 
linear quantities. For H 1 we find: 

H 1 « v\ + vf c + KY 1 (10.3.1) 

And: 

vjj(r) = | // c (r,r>i 1 (r',t)dV (10-3.2) 

£n/(r) = - J pHr,rU)u y (|r-r'|)/"(r')dV 
Plugging all this into Eq. (10.2.10): 
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ih^&r.i) * (j" (i^i +// c (r,r')}n 1 (r',t)d 3 r' 



+ v\r,t)^° k {r)- J pHr.r'.VuyQr -r'|)^(r')dV 
+ (/?[p°]-e k )^(r,t) 



(10.3.3) 



Plugging in the 1 st order definitions of n 1 and p 1 we find: 



ih — xplir.t) ~2Re) W k] (r,r')i}j(r',t)d 3 r' 




+ (H[p°] - e k )xpl(r, t) + v\r, t)^(r) 



(10.3.4) 



1 



^WjV) - uy(|r - r'|)^(r')^°(r) 



W kj (r,r') = 



_\r — r'\ 



+ flc(r,r')xp° k (r)xpf(r') 



We now take a crucial step to get rid of all the integrals and special 
dependence by expanding in terms of the orthonormal complete 
eigenfunction of the Hamiltonitan H[p ]. Because the evolution according to 
the Kohn-Sham equation preserves orthonormality. Thus, the response 
functions are orthogonal to all occupied MO's. This allows us to write: 

M 

iptir.t) = £ (4 s (t) + i4;(t))^ s °(r) (10.3.5) 

s=N + l 

We use the following index notations: 



j, k — occupied orbitals 



q,s — unoccupied orbitals 



p,n — all orbitals 



The W operator can be turned into a 4-index matrix: 
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(10.3.6) 

^(r)\p®(r')d 3 rd 3 r' 



W kpJn = jj W kj (r,r')^(r)^(r')d 3 rd 3 r' 

-u y (|r-r'|)^(r')^ (r) 
-> W kj {r,r') = £ W kpJn ^(r)K(r') 

pn 

Using these expansions we find: 
^m(4 s (t) + i4' s (t)> s °(r) 

s 

= 2 £ c; s (t)W kpJs Vp(r) (10.3.7) 

sJ,P 

s 

with: ftc*^ = e s — Now multiply by ipq(. r ) an d integrate over r and obtain 
a time-dependent algebraic equation: 

ih (c' kq {t) + ic^(t)) 

= 2 Z Z ^' s ^ w + ( c ^ (t) + ic ^^) ^ ( io - 3 - 8 ) 

Where: 

= / V q (T)v\r,t)r k (rWr (10.3.9) 
We can separate to two real equations by the real and imaginary parts: 
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-hc'; q {t) = lY^Wkq.jsCjsit) + C' kq (t)h0) qk + Vl q (t) 

s J (10.3.10) 

Cfcq(0 = C kq (t)a) qk 

This equation can be used to get a time-dependent response. It will be useful 
for describing the response to a complicated weak pulse v qk (t). 

C. Frequency-domain response 

Instead of viewing the response in time we can view it in frequency. This will 
allow us to address questions of response to sinusoidal fields, such as laser 
absorption. We Fourier- transform the coefficients: 

p CO 

d kq ((o) = c kq (t)e ia>t dt (10.3.11) 
Jo 

Note that c kq (0) = and so: — io)d kq (y) = f°° c kq (t)e l0)t dt. 

Note further that d' and d" are not real anymore. In fact, because c's are real 
we now have: d kq (oi)* — d kq (—oS). 

Applying a Fourier transform to the time-dependent equations gives: 

ihcod'^io)) = 2^ W kqJs d'j S (a)) + d' kq (co)hco qk + vl q {u>) 

i s (10.3.12) 

-i(x)d' kq {(x}) = d' k \{(D)(D qk 

Combining the equations we rid ourselves of d' and obtain an equation for 
d": 

2_^ 2h ~ 1(i) qk^kq,is + {^qk ~ ^ 2 )5{ qk ){ s j)]d kq {.(x)) = ih^cov^ (<o), (10.3.13) 

js 

The number of these coupled algebraic equations is infinite. However, in any 

practical application we use a finite number of states and so we have a finite 

number of equations. If M is the total number of states used we then for any 
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given perturbation one can get an approximate response by solving one of 
these two equations for dj s (a>). This involves solving a linear NM *NM set of 
equations and obtain the response to the perturbation. We can write it in 
matrix form: 

(A - co 2 0d"(co) = ih^oivKo)). (10.3.14) 

With: 

A q kjs = 2^ _1 a»q /c W /c£?JS + 0Jq fc 5( (J ifc)(sy) (10.3.15) 

We now prove that as long as a) qk for all possible q and k the matrix A has 
only real eigenvalues. In fact, a) qk is always positive unless there is a 
degeneracy at the HOMO and the level is not full. The eigenvalue equation is: 

(2h- 1 a) qk W kqJs + a) 2 qk 8 {qkXsj) )y S j = n 2 y qk (10.3.16) 

where y s j are the eigenvector components. Under the assumption we can 
bring this equation to the following equivalent form: 

^^(2h~ 1 W kq j s Ja) q ~ k ~a)^j + u q k8(qkXsj)) x sj = & 2x qk (10.3.17) 



JS 



Where x qk = -p=. These new equations involve the eigenvalues of a 

symmetric matrix. Thus, the eigenvalues H 2 are all real. Eq. (10.3.17) is called 
the Casida equations. 

D. Excitation energies from the algebraic 
treatment 

The algebraic equation for linear response derived in Eq. (10.3.14) allows 
calculation of a specific response problem. One incredible feature of this 
equation is that the matrix A on the left is independent of the perturbation v 1 . 
We can thus analyze A to get general response properties of the system, 
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without reference to any specific perturbation. We can define a "response 
matrix": 

R(co) = ih^oCA - a) 2 !)' 1 . (10.3.18) 
Where A is defined in (10.3.15)We can define right and left eigenvectors: 

Ay a = V-lVa, A T z a = a 2 a z a (10.3.19) 
It is easy to see from Eq. (10.3.16) that: 

*> qk (z a ) qk = (y a ) qk (10.3.20) 
We choose the normalizations so that: 

Zljfi = S a p = ^\za)js(y(i) js = ^ J ~ (ya)js(yp) js (10.3.21) 

is js SJ 

Collecting the vectors ina a matrix and the eigenvalues in a diagonal matrix 
H 2 we can write: 

A = y T a 2 z (10.3.22) 

And: 

RM =y T w^rf (10A23) 

We can now take to near an eigenvalue fl a . This has an infinite response. Thus 
we can then write: 

i/i _1 a> 

R(co) = y r a _^ z a (x) — > £l a (10.3.24) 

A general technique to study the resonance is to insert an infinitesimal 
imaginary part which we eventually take to zero: o> = £i a + irj. Then, after a 
few algebraic manipulations (i?(w) = y T a (n + Jx n } z a ~* R (.®-a + ir i) = 
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«(n. +a ,) = »- y l ^+g «. (10.3.25) 

Thus: 



lim W[R(£l a + irj)] = oo 

?j->oo 



T -1 

lim 3[i?(n a + ir])\ = y l a -rr—z a 
The imaginary part of the response is finite at resonances. 



(10.3.26) 
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